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SUMMARY 


Part  I:  A  Survey  of  Methods  for  Dynamic 
System  Identification 

This  section  of  the  report  summarizes  the  findings  of  the 
first  phase  of  a  continuing  study  concerned  with  the  utilization 
of  computers  in  the  determination  of  mathematical  models  for 
dynamic  systems.  In  this  initial  phase,  a  careful  review  of 
various  types  of  models  which  have  been  previously  used  or 
suggested  for  this  purpose  has  been  completed.  A  discussion 
of  these  models  and  the  associated  experimental  techniques  makes 
up  the  first  part  of  this  section.  Following  a  critical  evaluation 
of  the  earlier  methods,  a  new  formulation  of  the  general  model 
inference  problem  is  presented;  this  new  approach  has  been  named 
the  "parameter  space  method".  The  computational  advantages  of 
parameter  space  methods  are  discussed,  and  a  foundation  is  laid 
for  the  development  of  explicit  computational  procedures.  The 
continuation  of  this  research  will  include  the  investigation  of 
several  specific  techniques  for  achieving  the  "identification" 
of  nonlinear  systems. 
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Part  II:  Research  in  Optical  Coherence 

A  detailed  study  of  certain  aspects  of  the  theory  of  optical 
coherence  has  been  carried  out.  A  major  objective  of  the  study  is 
to  obtain  criteria  for  optimizing  optical  systems  employing  coherent 
sources  and  amplifiers  such  as  lasers.  This  is  the  first  in  a  series 
of  studies  to  investigate  the  effect  of  various  optical  system  devices 
such  as  antennas,  filters,  modulators  and  demodulators,  etc.,  on  the 
coherence  properties  of  optical  signals.  The  present  report  concerns 
itself  primarily  with  spatial  and  temporal  filters. 

A  general  coherence  function  is  defined  as  the  autocorrelation 
function  for  field  strengths.  The  field  strength  is  a  function  of 
three  space  variables  and  time.  The  coherence  function  is  therefore 
dependent  jointly  on  the  space  separation  and  the  time  separation 
of  the  measured  field  strengths.  It  is  also  a  function  of  the  origin 
in  space  and  time  if  the  process  is  nonstationary  in  those  variables. 
Limiting  cases  of  pure  spatial  coherence,  pure  temporal  coherence, 
pure  spatial  incoherence,  and  pure  temporal  incoherence  are  considered. 
The  concept  of  perfect  coherence  in  a  particular  independent  variable 
is  extended  to  include  all  deterministic  functions  which  are  completely 
specified  for  the  entire  range  of  that  variable. 

The  general  integral  equation  relating  image  and  object  field 
strengths  in  all  variables  is  defined.  A  technique  for  determining 


-iii- 


the  spatial  and  temporal  weighting  function  for  a  linear  optical 
transducer  is  suggested.  This  evaluation  is  complicated  by  the 
fact  that  the  diffraction  field  for  an  arbitrary,  aperture-field, 
distribution  has  only  been  found  (approximately)  as  the  steady-state 
response  to  a  monochromatic  source.  Since  temporally  incoherent  or 
partially  coherent  sources  are  common,  the  weighting  function  in 
both  space  and  time  must  be  found.  This  function  may  be  found 
approximately  from  the  steady-state,  sinusoidal,  response  by  taking 
the  inverse  Fourier  transform  of  the  spatial  response  to  a  sinusoid 
over  the  frequency  range  for  which  the  sinusoidal  response  function 
is  valid.  This  inversion  is  complicated  by  the  fact  that  the  spatial 
response  is  a  function  of  the  temporal  frequency. 

The  technique  of  obtaining  the  spatial  and  temporal  weighting 
function  makes  possible  the  evaluation  of  both  the  transient  and 
steady-state,  statistical,  behavior  of  far-field  diffraction  patterns 
for  sources  exhibiting  partial  coherence  in  both  space  and  time.  In 
particular,  information  rates  may  be  calculated,  and  the  transducer 
may  be  optimized  with  respect  to  its  variable  parameters  by  maximizing 
the  image  information  rate  with  respect  to  these  parameters. 
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1.  INTRODUCTION 


1.1  Purpose  and  Scope 

This  report  summarizes  the  findings  of  the  first  phase  of  a  continuing 
research  program  concerned  with  the  utilization  of  computers  in  the 
determination  of  mathematical  models  for  dynamic  systems  from  experimental 
observations.  In  this  initial  phase,  a  careful  review  of  various  types 
of  models  which  have  been  previously  used  or  suggested  for  this  purpose 
has  been  completed.  A  discussion  of  these  models  and  the  associated 
experimental  techniques  makes  up  the  first  part  of  this  report.  Following 
a  critical  evaluation  of  the  earlier  methods,  a  new  formulation  of  the 
general  model  inference  problem  is  presented;  this  new  approach  has  been 
named  the  "parameter  space  method".  It  is  believed  that  parameter  space 
techniques  offer  significant  advantages  for  the  actual  computational 
determination  of  models  for  nonlinear  systems  in  practical  situations. 

The  continuation  of  the  research  reported  herein  is  aimed  at  the  develop¬ 
ment  of  explicit  computational  algorithms  for  achieving  the  identification 
of  nonlinear  systems  by  parameter  space  operations.  This  research  will 
provide  the  basis  for  future  reports  dealing  with  this  topic. 

1.2  Identification  of  Dynamic  Systems 

’  Since  the  time  of  Newton  it  has  been  known  that  the  dynamical 
behavior  of  mechanical  systems  is  governed  by  differential  equations. 

In  the  ensuing  centuries  since  Newton's  discoveries,  a  complete  theory 
of  classical  mechanics  has  been  constructed  which,  in  principle,  permits 
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the  description  of  the  motion  of  any  rigid  body  system  in  terms  of  a  set 
of  coupled  differential  equations.  In  a  similar  fashion,  electric  circuit 
theory  has  provided  differential  equation  descriptions  for  lumped  constant 
electrical  devices.  More  generally,  it  is  now  common  knowledge  that  a 
great  many  processes  involving  storage  and  interchange  of  mechanical, 
electrical,  and  other  forms  of  energy  are  properly  described  by  complicated 
sets  of  differential  equations  (or  partial  differential  equations  if  the 
systems  are  distributed  rather  than  lumped).  Systems  of  this  type  are 
usually  referred  to  as  "dynamic"  systems  in  engineering  literature  (1)\ 

While  basic  physical  theories  often  permit  the  form  of  the  differen¬ 
tial  equation  for  a  dynamic  system  to  be  written  down,  as  a  rule  the 
parameters  of  a  particular  system  can  be  determined  only  by  measurement. 

For  example,  the  equation  for  a  pendulum  with  viscous  damping  can  be 
obtained  from  the  simple  moment  equation: 

10  +  BO  =  -  M  gl  sin  0  (l) 

However,  before  any  analysis  of  the  behavior  of  a  particular  pendulum  may 
be  effected,  it  is  necessary  to  determine  the  constants  in  the  equation. 

Such  determination  of  a  specific  quantitative  model  has  been  variously 
referred  to  as  the  "identification  problem"  (2),  the  "characterization 
problem"  (3),  and  the  "parameter  estimation"  problem  (4), depending  somewhat 
upon  the  methods  employed;  it  can  be  accomplished  only  by  experimentation. 


In  this  report,  superscript  numbers  are  used  to  indicate  footnotes  while 
bracketed  numbers  refer  to  the  list  of  references  at  the  end  of  the  text. 
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In  the  case  of  a  device  as  simple  as  a  pendulum,  the  evaluation  of  the 
unknown  constants  is  not  difficult  unless  a  high  degree  of  precision  is 
required.  On  the  other  hand,  there  are  many  practical  situations  in  which 
the  measurement  of  system  constants  is  extremely  difficult  due  either  to 
the  complexity  of  the  processes  involved  or  the  subtleness  of  the  effects 
to  be  measured.  Furthermore,  even  after  the  system  parameters  have  been 
established  to  a  sufficient  accuracy,  analytic  solution  of  the  resulting 
family  of  differential  equations  is  an  extremely  difficult  task  in  all 
but  the  simplest  of  situations. 

For  the  reasons  just  mentioned,  engineers  and  physicists  over  the 
past  several  decades  have  sought  ways  of  characterizing  physical  systems 
other  than  by  a  set  of  differential  equations.  Such  efforts  have  been 
quite  successful  with  respect  to  the  important  but  restrictive  class  of 
systems  possessing  the  properties  of  linearity  and  time  invariance.  For 
more  general  classes  of  systems,  results  have  been  sparse.  If  the  require¬ 
ment  of  time  invariance  i6  relaxed  while  linearity  is  retained,  weighting 
functions  offer  an  alternative  to  differential  equation  descriptions  although 
they  may  be  very  difficult  to  obtain.  For  general  nonlinear  systems,  it  is 
only  by  considering  a  system  as  an  operator  effecting  a  transformation  from 
an  input  function  space  to  an  output  space  that  it  has  been  possible  to 
arrive  at  alternate  descriptions. 

The  following  paragraphs  discuss  various  methods  which  have  been  used 
to  describe  dynamic  systems  and  suggest  a  particular  characterization  to  be 
used  as  the  basis  for  the  development  of  practical  computation  procedures 
for  identifying  unknown  nonlinear  systems.  The  discussion  will  be  confined 
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to  lumped  constant  systems;  i.e.,  to  systems  described  by  ordinary  rather 
than  partial  differential  equations.  For  the  sake  of  completeness  and 
continuity!  the  development  to  follow  begins  with  a  review  of  methods  which 
have  been  developed  exclusively  for  the  characterization  of  linear  systems. 

2.  LINEAR  SYSTEM  IDENTIFICATION 

2.1  Frequency  Response 

2.1.1  Introduction 

When  a  system  is  linear,  the  "principle  of  superposition"  may  be 
applied  to  decompose  arbitrary  forcing  functions  and  initial  conditions 
into  component  parts  whose  effects  may  be  more  easily  analyzed.  Hie 
individual  responses  of  these  components  may  then  be  added  to  produce  the 
total  system  response.  If  the  system  under  consideration  is  time  invariant 
as  well  as  linear,  this  type  of  analysis  is  especially  simple  and  effective. 

In  particular,  when  a  sinusoidal  or  complex  exponential  decomposition  of 
the  input  signal  is  utilized,  Fourier  transform  techniques  may  be  applied 

to  determine  the  system  output.  In  this  approach,  the  system  under  considera- 

2 

tion  is  completely  determined  by  the  function,  Q(«),  defined  by 

Y(w)  =  G(«)  X(w)  (2) 

where  X(ai)  and  Y(a>)  are  Fourier  transforms  of  the  input  function,  x(t), 
and  the  output  function,  y(t),  respectively.  The  function  G(«)  is 
usually  called  the  system  "frequency  response"  (5,6,7).  For  a  linear  time 

2 

While  input  and  output  variables  will  be  treated  as  scalars  in  this 
discussion,  no  difficulties  are  experienced  if  equation  1  is  a  vector- 
matrix  relationship. 
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invariant  system,  x(t)  and  y(t)  are  related  by 


Ay(t)  *  Bx(t)  (3) 

where  A  and  B  are  linear  constant  coefficient  differential  operators; 
the  frequency  response,  G(u),  is  therefore  a  rational  polynomial  in  o> 
in  such  circumstances. 

2.1.2  Experimental  Determination  of  Frequency  Response 

Within  the  class  of  stable,  linear,  time  invariant  systems,  experimen¬ 
tal  evaluation  of  the  frequency  response  function  is  carried  out  by  a 
procedure  which  is  (mathematically)  independent  of  the  nature  of  the  system. 
This  is  a  major  advantage  of  the  frequency  method.  No  assumptions  regarding 
the  system  must  be  made  beyond  the  basic  constraint  that  it  be  stable, 
linear,  and  time  invariant.  In  the  Usual  technique  employed  to  experimen¬ 
tally  determine  the  frequency  response,  G(a>),the  input,  x,  is  chosen  to 
be  the  function, 


x(t)  =  a  cos  w^t  (4) 

so  that  after  a  time  sufficient  for  the  decay  of  transients  has  elapsed 

y(t)  *  b  cos  (u>^t  +  0) .  (5) 

It  is  easily  demonstrated  that  the  "frequency  response"  at  o>  =  u>^  is  then 
given  by 
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a  +  JP 


(6) 


gC^) 


b_ 

a 


The  variables  ■—  and  0  are  commonly  called  "gain"  and  "phase  shift" 
respectively  and  are  measurable  by  standard  electronic  devices. 

The  "Nyquist  plot"  is  the  locus  of  Q(w)  in  a,P  coordinates. 
Special  purpose  devices  have  been  constructed  to  plot  Nyquist  diagrams 
automatically  since  the  point  by  point  determination  of  G((o)  may  be 
very  tedious  (8). 

2.1.3  Measurement  of  Frequency  Response  in  the  Presence  of  Interfering 
Noise 

When  extremely  precise  measurement  is  attempted  or  uncontrollable 
sources  of  interfering  noise  are  present,  it  is  found  that  the  measured 
gain  and  phase  shift  are  random  variables.  To  circumvent  this  difficulty, 
cross-correlation  techniques  may  be  used  to  discriminate  against  the  noise 
(9).  From  equation  6,  y(t)  may  be  written 


y(t)  «  aa  cos  u^t  -  aP  sin  u^t 


(7) 


With  noise,  n(t)  added,  this  becomes 


1 

y  (t)  *  aa  cos  ai^t  -  aP  sin  w^t  +  n(t) 


(8) 


Now  assuming  that  y(t)  and  n(t)  are  uncorrelated,  it  follows  that 

.T 


2  r  » 

lim  ^  J  y  cos  <»^t  dt  *  aa 


T  -  CD 

j  ( 

lim  |  J  -  y  sin  w^t  dt  «  aP 
T  *oo  C 


(9) 
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Finite  time  approximations  to  these  cross-correlation  functions  can  yield 
arbitrarily  precise  measurements  of  a  and  P  in  the  presence  of  arbitrarily 
large  noise  voltages. 

2.1.4  Extension  of  Frequency  Methods  to  Nonlinear  Systems 

As  a  result  of  the  simplicity  and  practical  utility  of  the  frequency 
response  characterization  of  linear  systems,  attempts  have  been  made  to 
extend  frequency  methods  to  certain  types  of  nonlinear  systems.  Perhaps, 
the  most  significant  of  these  is  the  "describing  function"  method  due  to 
R.  J.  Kochenberger  (10).  Kochenberger's  method  consists  of  determining 
the  gain  and  phase  characteristics  of  a  nonlinear  element  with  respect  to 
the  fundamental  component  of  the  output  only.  In  general,  the  frequency 
response  defined  in  this  way  depends  not  only  upon  the  frequency  of  the 
test  signal,  but  also  upon  its  magnitude.  Describing  functions  have  found 
considerable  use  in  the  investigation  of  limit  cycles  in  nonlinear  systems. 
Since  the  method  is  approximate,  it  is  of  relatively  little  value  in 
determining  the  transient  behavior  of  nonlinear  systems. 

The  method  for  experimentally  determining  describing  functions  is 
the  same  as  for  experimentally  determining  frequency  response  except  that 
the  describing  function  depends  on  both  a  and  ui  in  equation  6  rather 
than  on  oi  alone. 

2.2  Determination  of  the  System  Transfer  Function 

2.2.1  Relationship  Between  the  System  Laplace  Transform  and  the  Frequency 
Response  Function 

If  the  Laplace  transform  is  applied  to  both  sides  of  equation  3, 
the  result  is 

A(s)  Y(s)  .  B(s)  X(s)  (10) 
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so 

frfr  ■  xffr  -  H(s)  (11) 

The  function  H(s)  is  commonly  called  the  system  transfer  function. 

Again,  since  A  and  B  are  assumed  to  be  linear  constant  coefficient 
differential  operators,  H(s)  is  a  rational  polynomial.  As  is  well  known, 
the  frequency  response  function,  G(u),  may  be  obtained  from  H(s)  by 
utilizing  the  relation 

G(w)  -  H( jm)  '  (12) 

Furthermore,  since  H(a)  is  an  analytic  function,  its  behavior  along  the 
imaginary  axis  in  the  s  plane  completely  specifies  its  behavior  over  the 
entire  complex  plane.  So,  conversely,  experimentally  evaluated  frequency 
response  may  be  used  to  determine  H(s)  for  stable  systems.  The  process 
of  obtaining  H(s)  from  G(w)  is  called  the  "approximation  problem"  in 
network  theory  (11,12,13,14). 

2.2.2  Direct  Measurement  of  the  System  Transfer  function 

Up  to  the  present  time,  the  frequency  response  method  has  been  over¬ 
whelmingly  favored  in  practice  for  the  experimental  characterization  of 
linear  time  invariant  systems.  Recently,  however,  a  few  investigators 
have  devised  methods  for  measuring  pole  and  zero  locations  directly  by 
using  special  test  signals.  These  are  sometimes  referred  to  as  "time 
domain"  methods  (15,16,17,18).  While  these  techniques  have  the  advantage 
that  the  approximation  problem  is  by-passed  entirely,  the  test  equipment 
required  is  generally,  more  complicated.  In  any  event,  the  fined  result 
of  the  measurement  is  mathematically  equivalent  to  the  frequency  response. 
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2.3  Weighting  Function  Methods 

2.3.1  Relation  to  the  System  Transfer  Function 

The  behavior  of  any  linear  system  with  respect  to  an  input  or 
disturbance  at  a  particular  point  is  completely  determined  by  its  weighting 
function  or  "impulse  response"  (19).  If  h(t,T)  is  the  weighting  function 
associated  with  a  linear  differential  operator,  L,  and  y(t)  are  input 
and  output  variables  respectively,  then  the  solution  to 

Ly  -  x  (13) 

is  given  by 

-  t 

y(t)  -  I  h(t,r)  x(t)  dr 
•  00 

If  L  is  also  time  invariant,  then  h(t,T)  becomes  a  function  of 
only  and  is  related  to  the  system  transfer  function  by 

{h(t-T)}  -  e~TS  H(s)  (15) 

where  denotes  the  Laplace  transform  operator  (19).  Thus  the  weighting 
function  and  frequency  response  form  a  Laplac e-Fourier  transform  pair. 

While  the  weighting  function  bears  a  simple  relationship  to  frequency 
response  only  in  the  time  invariant  case,  it  is  nevertheless  (in  light  of 
equation  14)  a  perfectly  general  way  of  characterizing  any  linear  system. 

The  determination  of  h(t,r)  from  the  operator  L  is  apt  to  be  a  formidable 
task,  however.  Most  often,  numerical  methods  or  analog  simulation  are 
required  to  find  h(t,t). 


(14) 

(t-r) 
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2.3*2  Direct  Measurement  of  Weighting  Functions 

If  the  system  under  test  is  initially  in  a  relaxed  state,  then  for 
t  >  0,  the  lower  limit  in  the  integral  of  equation  14  may  be  replaced  by 
zero.  If  x(t)  is  then  chosen  to  be  a  narrow  pulse  of  unit  area  centered 
at  t  -  t^  >  0,  the  system  weighting  function  can  be  measured  directly. 
Suppose  that  x(t)  is  narrow  enough  to  be  sufficiently  approximated  by  a 
unit  impulse  function,  6(t-t^).  Then 

y(t)  ■  J  h(t,T)  j(t-t1)  dr  -  hU,^)  (16) 

By  repeating  this  experiment  with  various  values  of  ^  ,  a  family  of  curves 
for  h(t,r)  can  be  constructed. 

When  L  is  time  invariant,  it  is  sufficient  to  make  a  single  exper¬ 
iment  with  t^  *  0  since  h(t,r)  depends  only  on  t-T.  In  this  respect, 
weighting  function  measurement  of  system  dynamics  is  markedly  superior  to 
the  frequency  method  which  is  a  point  by  point  procedure.  On  the  other 
hand,  the  equipment  needed  to  record  transient  behavior  is  usually  more 
complicated  and  restricted  in  applicability  than  that  used  in  determining 
frequency  response. 

E.  Miskin  and  R.  A.  Haddad  have  described  an  adaptive  system  which 
uses  integrals  of  impulse  functions  (i.e,,  steps,  ramps,  and  parabolas) 
for  both  identification  and  control  of  an  object  whose  weighting  function 
is  initially  unknown  (20).  This  method  permits  the  control  and  identifica¬ 
tion  functions  to  proceed  simultaneously  but  it  is  restricted  to  situations 
where  measurement  noise  is  insignificant. 
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2.3.3  Determination  Of  Weighting  Functions  by  Cross-Correlation 

All  of  the  methods  mentioned  thus  far  for  characterizing  linear 
systems  share  a  common  weakness.  In  order  to  make  measurements  on  the 
device  in  question,  it  is  necessary  to  remove  it  from  its  normal  function 
and  apply  special  test  signals.  That  is,  the  methods  described  so  far 
involve  essentially  laboratory  or  "bench  test"  types  of  measurements. 

There  are  many  important  situations  in  which  it  is  either  impossible  or 
undesirable  to  create  such  controlled  conditions.  In  such  cases,  cross- 
correlation  of  input  and  output  may  be  used  to  discriminate  against  the 
unwanted  effects. 

The  application  of  croBS-correlation  to  frequency  testing  has  already 
been  discussed  in  this  report.  In  1950,  J.B.  Wiesner  and  Y.  W.  Lee  pointed 
out  that  cross-correlation  can  also  be  used  to  advantage  in  the  measurement 
of  weighting  functions  (21).  Following  their  suggestion,  suppose  that  the 
input  to  a  system,  x(t),  is  a  wide-band  stationary  Gaussian  process.  Then 
the  autocorrelation  function  of  this  process  is  given  by 

j?  (t)  ■  E  (x(t)  x(t+T ) }  *  N  6(t)  ■  0  (-r)  (17) 

XX  O  XX 

where  Nq  is  the  (two  sided)  spectral  density  of  the  input  variable,  x(t)  . 
If  the  cross-correlation  between  the  input,  x,  and  the  output,  y,  is 
defined  by 


ff  (t,r)  »  E  {y(t)  x( t+T ) ) 


(18) 
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then 


0yx(t*-T)  -  E  (y(t)  x(t-r)} 


=  E  [x(t-T)  h(t,a)  x(a)  da} 
Jo 

Jt 

=  |  E  (x(t-r)  x(a)}  h(t,a)  da 
Jo 
/•t 

«  j  Nq  4(a-t+r)  h(t,a)  da 


N  h(t,t-t) 
o 


(19) 


or 


h(t,t-r)  *  — 


N 


In  the  special  situation  when  0  is  not  time  dependent,  then 

yx 


(20) 


equation  20  reduces  to 


0  C-t) 

h(r)  »  h(--r)  »  - 


(21) 


which  may  be  replaced  by  a  time  average: 

i  rT 


h( 


lim  -■  F  y(t)  x(t-T)  dt 
t)  -  T  Jo 


(22) 


T  a: 


This  is  the  result  originally  pointed  out  by  Lee  and  Weisner  (21)  and  the 
one  most  often  used  in  weighting  function  methods  for  system  identification. 
The  input  process,  x(t),  may  represent  either  normal  operating  signals  alone 
or  such  signals  plus  a  deliberately  added  wideband  test  signal.  In  the  event 
that  a  special  test  signal  is  added,  it  is  only  this  signal  which  ia  cross- 
correlated  with  the  system  outjwt. 
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2.3.4  Experimental  Evaluation  of  Cross-Correlation  Functions  for  Time 

Invariant  Linear  Systems 

Following  the  original  publication  of  Lee  and  Weisner,  a  sizable 

number  of  investigators  have  explored  various  techniques  for  the  experimental 

evaluation  of  0  (t).  Despite  the  considerable  originality  of  some  of 

yx 

these  methods,  nearly  all  of  them  contain  the  basic  elements  of  delay, 
multiplication,  and  averaging  as  shown  in  Figure  1.  Since  the  averaging 
time,  T,  must  be  finite,  the  output  of  the  process  identifier  is  h(x), 
an  estimator  of  h(t).  The  mean  square  deviation  of  h(-r)  from  h(-r)  may 
be  made  as  small  as  desired  by  increasing  the  averaging  time  sufficiently. 

The  next  few  paragraphs  discuss  some  specific  hardware  mechanizations  of 
equation  22  as  well  as  some  minor  variations  on  the  basic  scheme. 

J.  A.  Aseltine,  et  al.  have  suggested  that  a  binary  test  signal  offers 
particular  advantages  in  the  experimental  determination  of  weighting  func¬ 
tions  (22).  If  an  input  is  used  which  switches  from  one  level  to  the  other 
in  the  manner  of  a  Poisson  process,  then  the  delay  time,  T,  is  easily 
realized  by  either  a  delay  line  or  a  shift  register  while  the  multiplication 
required  is  replaced  by  a  simple  on-off  electronic  gate.  While  a  Gaussian 
random  process  has  most  often  been  used  as  a  probing  signal  for  weighting 
function  determination,  equation  21  remains  valid  for  any  input  with  an 
autocorrelation  function  which  is  sharp  compared  to  h(-r).  Consequently, 
a  binary  signal  is  entirely  satisfactory  as  a  test  signal  so  long  as  the 
average  number  of  switchings  per  second  considerably  exceeds  the  bandwidth 
of  the  system  under  test. 


-13- 


Q.  L.  Turin  has  discussed  the  replacement  of  explicit  elements  of 
Figure  1  by  a  filter  matched  to  a  particular  test  signal  (23).  When  this 
is  done,  the  physical  equipment  involved  is  greatly  simplified.  Moreover, 
h(t)  is  obtained  as  a  continuous  function  of  time  rather  than  on  a  point 
by  point  basis  as  in  the  conventional  cross-correlation  method.  On  the  other 
hand,  Turin's  method  does  not  permit  long  averaging  times  to  be  used  to 
discriminate  against  outputs  resulting  from  extraneous  noise  or  normal 
operating  signals.  W.  W.  Lichtenberger  has  suggested  that  this  difficulty 
can  be  circumvented  by  averaging  the  filter  response  over  a  number  of 
measurements  (24),  However,  this  again  yields  a  point-wise  estimate  of 
h(T)  unless  an  analog  memory  is  available  to  store  successive  samples  of 
h(t)  in  their  entirety. 

Several  investigators  have  used  multiple  sinusoidal  test  signals 
followed  by  syncronous  detection  in  place  of  a  Gaussian  input  (25,26,27,28). 

In  these  schemes  it  is  system  coefficients  which  are  to  be  determined  directly 
rather  than  the  weighting  function.  Unlike  simple  frequency  response  testing, 
these  methods  use  averaging  following  the  syncronous  detection  and  are  there¬ 
fore  able  to  function  in  the  presence  of  interfering  signals.  Despite  the 
fact  that  the  test  signals  are  sinusoidal,  this  type  of  measurement  is  a 
correlation  method  which  fits  into  the  framework  of  Figure  1. 

Quite  a  number  of  investigations  have  been  directed  toward  the 
problem  of  determining  weighting  functions  by  correlation  methods  utilizing 
normal  operating  records  without  the  introduction  of  special  test  signals. 

The  first  significant  work  in  this  direction  appears  to  have  been  accomplished 
by  T.  P,  Goodman  and  J.  B.  Seswick  (29).  They  constructed  a  special  piece 
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of  equipment  involving  a  tapped  delay  line  and  manually  set  weighting 
potentiometers  to  match  the  cross-correlation  properties  of  a  system  under 
test.  The  unknown  weighting  function  evaluated  at  evenly  spaced  time 
intervals  appeared  as  potentiometer  settings  at  the  end  of  an  iterative 
manual  adjustment  process.  R.  £.  Kalman  has  demonstrated  an  automatic 
iterative  computational  technique  for  obtaining  a  least  squares  estimate 
of  a  z- transformed  version  of  the  system  weighting  function  during  normal 
operation  (30).  P.  Joseph  et  al.  have  proposed  a  variation  in  Kalman's 
method  which  allows  initial  conditions  to  be  considered  (31).  R.  B.  Kerr 
and  W.  H.  Surber  have  examined  the  relationship  between  time  of  observation 
and  accuracy  of  weighting  function  estimation  (32).  V.  V.  Solodovnikov 
and  A.  S.  Uskov  have  discussed  transform  methods  for  the  solution  of  the 
equation 

*yx(T)  =  J  **xx(T-t}  h(t)  dt  '  (23) 

for  h(t)  (33).  This  equation  arises  when  the  input  correlation  function  is 
not  extremely  narrow  compared  to  h(t)  and  the  averaging  time  of  Figure  1 
is  allowed  to  become  infinite. 

There  have  been  many  other  publications  dealing  with  cross-correlation 
since  the  original  observation  of  Lee  and  Weisner.  Since  the  present  discus¬ 
sion  is  intended  simply  to  place  the  methods  proposed  in  this  report  in 
perspective!  the  examination  of  cross-correlation  techniques  has  not  been 
exhaustive.  The  bibliographies  attached  to  the  publications  referenced 
provide  many  other  sources  of  information  concerning  this  topic. 
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2.5*5  Identification  of  Tine  Varying  Linear  Systems  by  Cross-Correlation 

The  experimental  determination  of  time  varying  weighting  functions 

using  equation  20  is  seriously  hampered  by  the  fact  that  expected  value 

operators  cannot  be  replaced  by  time  averages.  Rather,  it  is  necessary  to 

estimate  0  (t-r)  by  averaging  the  results  of  many  separate  experiments, 

yx 

Furthermore,  since  h(t,*r  )  is  a  function  of  two  variables,  these  repeated 
experiments  must  be  carried  out  for  each  value  of  t  to  be  considered. 

Thus  the  total  number  of  experiments  involved  is;  very  large.  Moreover, 
the  system  must  be  restored  to  a  reference  condition  (t«0)  at  the  start 
of  each  experiment.,  It  may  not  be  feasible  to  achieve  this  type  of  operation 
in  many  systems  of  practical  concern.  Despite  these  difficulties,  equation 
20  appears  to  offer  promise  of  useful  application  in  special  situations. 

The  necessary  storage  and  averaging  could  be  accomplished  by  a  digital 
machine  possessing  analog  inputs  and  outputs.  As  far  as  is  known  to  the 
writer,  no  application  of  this  approach  has  been  made  to  date. 

2.3.6  Impulse  Response  Measurement  by  Regression  Analysis 

Due  to  the  force  of  tradition  more  than  anything  else,  analog 
measuring  and  computing  devices  have  been  assumed  almost  exclusively  in  the 
investigation  of  practical  procedures  for  impulse  response  measurement. 
Mathematically,  this  means  that  the  techniques  used  have  been  restricted  to 
processes  of  analysis  to  the  virtual  exclusion  of  algebraic  methods.  By 
contrast,  classical  statistics  is  concerned  mainly  with  discrete  data 
obtained  from  repeated  experiments  so  algebraic  methods  predominate  in 
that  subject.  Thus,  when  conventional  statistical  methods  are  applied  to 
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the  measurement  of  a  system  impulse  response)  a  fresh  viewpoint  results. 

An  example  of  the  statistical  approach  has  been  provided  by  M.  J, Levin 
in  an  investigation  dealing  with  the  estimation  of  impulse  response  at 
discrete  values  of  time  by  linear  regression  analysis  (34).  In  Levin's 
approach,  the  convolution  integral 

y(t)  -  F  h(r)  x(t-r)  dT  (24) 

Jo 

is  replaced  by  an  approximating  convolution  summation 

N 

y(NT)  -  £  h(nT)  {  x  [(N-n)T]}  T  (25) 

n-0 

In  this  expression,  x  and  y  are  physically  measurable  while  h(nT)  is 
unknown.  Since  the  equation  is  linear  in  the  unknown  weighting  function 
values,  a  set  of  N+l  simultaneous  equations  resulting  from  application 
of  equation  25  may  be  solved  for  h(nT),  n»0,  1,  N,  by  matrix  inversion. 
However,  since  measurement  errors  are  invariably  present,  Levin  suggests 
that  redundant  data  be  taken  and  h(nT)  be  determined  by  least  squares 
regression  analysis.  It  is  shown  in  his  paper  that  this  results  in  estimates 
which  are  optimum  in  several  senses. 

By  substituting  equation  25  for  equation  24  the  difficult  problem  of 
finding  an  Inverse  operator  is  reduced  to  the  much  easier  one  of  finding  the 
inverse  of  a  matrix.  This  is  an  operation  well  suited  to  a  digital  computer. 
In  a  paper  dealing  with  future  trends  in  engineering  analysis,  Denis  Gabor 
has  pointed  out  that  it  is  quite  typical  that  the  discrete  problem  should 
yield  a  solution  more  readily  than  the  continuous  problem  (35) •  It  is 
further  stated  in  his  article  that  algebraic  methods  can  be  expected  to 
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replace  analytic  approaches  to  a  great  many  problems  in  applied  mathematics. 
This  opinion,  it  seems,  is  supported  by  the  fact  that  high  speed  computers 
are  now  available  to  most  mathematicians  and  engineers.  While  these  machines 
have  great  algebraic  power,  they  are  not  naturally  suited  to  limiting  processes. 

The  point  of  view  expressed  by  Levin  appears  to  offer  an  attractive 

alternative  to  correlation  methods.  The  continuation  of  the  research  die- 
* 

cussed  in  this  report  involves  the  application  of  similar  approaches  to 
nonlinear  systems. 

3.  FUNCTION  SPACE  DESCRIPTIONS  FOR  DYNAMIC  SYSTEMS 
3.1  Introduction 

All  of  the  methods  employed  for  linear  system  identification  appear 
to  run  aground  when  attempts  are  made  to  apply  them  to  general  nonlinear 
systems.  This  unfortunate  situation  results  from  the  necessity  of  abandoning 
the  principle  of  superposition  in  dealing  with  nonlinear  systems.  In  order 
to  gain  a  vantage  point  on  nonlinear  systems  comparable  to  the  conventional 
treatment  of  linear  systems,  it  seems  to  be  necessary  to  introduce  the  much 
more  sophisticated  idea  of  transformations  defined  over  function  spaces  (36). 
That  is,  it  is  necessary  to  recognize  that  the  output  of  a  physical  system 
in  the  most  general  case  depends  upon  the  entire  past  history  of  its  input 
in  some  complicated  nonlinear  fashion.  This  section  provides  a  brief  summary 
of  the  function  space  point  of  view. 

The  application  of  function  space  methods  is  considerably  facilitated 
by  the  expansion  of  input  functions  into  orthogonal  series.  The  particular 
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type  of  orthogonal  functions  used  for  this  expansion  will  depend  upon  the 
nature  of  the  input  and  the  test  being  conducted.  For  example,  if  the 
input  is  periodic!  then  a  Fourier  series  expansion  is  appropriate.  On 
the  other  hand,  if  the  input  is  not  known  to  be  periodic  but  is  known  to 
have  finite  bandwidth,  B,  it  can  be  completely  represented  by  tine  samples 
separated  by  intervals  equal  to  1/2B.  This  representation  is  called  a 
"signal  space"  expansion  in  communication  theory  (37).  If  the  input  is 

a  general  random  process  unspecified  except  for  the  properties  of  contin- 

/ 

uity  and  finite  variance,  then  the  entire  past  history  of  its  values  may 
be-  summarized  by  the  coefficients  of  an  expansion  in  terms  of  Laguerre 
functions  (38).  Whatever  the  type  of  expansion  chosen,  the  result  is  that 
the  input  may  be  thought  of  as  either  a  point  or  a  curve  in  a  space  of 
infinite  dimension;  i.e.  in  a  Hilbert  space.  Representation  as  a  point 
occurs  when  the  entire  input  is  known  in  advance  while  a  space  curve 
results  when  the  input  is  a  random  process  with  only  past  values  known 
and  a  "backward  looking"  expansion  is  used. 

When  the  input  to  a  system  has  been  appropriately  expanded,  the  out¬ 
put  of  the  system  may  be  regarded  as  the  result  of  a  mapping  or  transformation 
0  defined  over  the  input  space.  Physically,  this  transformation  is  produced 
by  the  reaction  of  the  system  under  test  to  a  "probing"  signal  applied  at 
an  inpiit  (38).  When  the  entire  input  can  be  represented  as  a  stationary 
point  in  a  Hilbert  space,  the  output  can  be  thought  of  as  another  point  in 
a  similar  space.  This  is  the  approach  used,  for  example,  in  steady  state 
harmonic  analysis  of  linear  systems.  When  the  input  is  represented  by  a 
moving  point  in  a  Hilbert  space  as  in  the  case  of  a  random  input  process, 


then  the  output  space  is  just  the  real  line.  That  is,  the  entire  past 
history  of  the  input  determines  the  present  output;  future  inputs  and  outputs 
are  not  known. 

When  the  function  space  point  of  view  is  adopted,  the  system  identifica¬ 
tion  problem  becomes  a  matter  of  determining  what  transformation  characterizes 
the  system  under  test.  Again,  there  is  considerable  freedom  of  choice  in 
selecting  a  method  for  describing  this  transformation.  For  example,  if  the 
system  is  known  to  be  linear  and  time  invariant,  then  the  frequency  response 
is  an  appropriate  characterization  of  the  system  input  to  output  transforma¬ 
tion.  Specifically,  if  an  input  x,  is  periodic  in  an  interval  T  and  is 
square  integrable  within  that  interval,  then 
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The  complex  function,  G(u>),  is  just  the  frequency  response  defined  by 
equation  2.  Since  x  may  be  thought  of  as  a  point  or  vector  in  a  space 
whose  coordinates  are  the  complex  quantities,  an,  and  y  is  likewise  a 
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point  in  the  bQ  space,  equation  28  is  an  example  of  a  particularly  simple 
transformation  from  a  Hilbert  space  to  a  Hilbert  space. 

While  equation  28  represents  merely  a  rephrasing  of  the  significance 
of  frequency  response,  as  indicated,  the  function  space  point  of  view  is  not 
limited  to  linear  systems.  It  is  possible  to  choose  families  of  orthogonal 
functions  defined  over  the  input  Hilbert  space  coordinates  which  are  capable 
of  representing  arbitrary  nonlinear  operators  ^ .  Norbert  Wiener,  for  example, 
chose  to  represent  a  Gaussian  input  process  in  terms  of  normalized  Laguerre 
functions  and  then  described  the  system  output  as  an  expansion  involving 
products  of  normalized  Hermite  functions  whose  arguments  are  the  Laguerre 
coefficients  (3,38).  In  electing  to  use  these  particular  series,  Wiener 
restricted  himself  to  the  use  of  broadband  Gaussian  random  processes  as  a 
source  of  probing  or  test  signals  and  to  time  invariant  operators  having 
the  property  that  inputs  applied  arbitrarily  far  in  the  past  have  an 
arbitrarily  small  influence  on  the  present  system  output.  This  latter 
restriction  means,  for  one  thing,  that  the  Wiener  expansion  cannot  be  used 
for  the  study  of  unstable  or  oscillating  nonlinear  systems. 

It  is  certainly  possible  to  choose  expansions  different  from  the  one 
chosen  by  Wiener  to  characterize  a  nonlinear  operator.  For  example,  A.  Bose 
has  discussed  an  expansion  which  involves  basically  the  idea  of  partitioning 
a  Hilbert  space  into  cells  and  associating  a  particular  output  with  each 
cell  (3).  A.  W.  Balakrishnan  has  discussed  the  application  of  polynomials 
defined  over  Hilbert  spaces  to  the  problem  of  nonlinear  operator  representa¬ 
tion  (39).  L.  A.  Zadeh  has  provided  a  tutorial  exposition  of  the  Hilbert  . 

^  The  operator  may  be  of  an  arbitrary  nature  providing  only  that  it 
produces  outputs  which  lie  in  the  selected  output  function  space. 
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space  characterization  of  nonlinear  operators  wherein  several  more  expansions 
are  described  (36) 

No  matter  what  orthogonal  functions  are  used  as  a  basis  for  the  representa¬ 
tion  of  the  input  and  the  system  transformation,  with  each  such  pair  of 
expansions  there  will  be  a  class  of  inputs  and  a  class  of  transformations 
which  can  be  represented.  The  choice  of  representation  is  thus  determined 
by  the  problem  to  be  solved. 

3.2  A  Function  Space  Definition  of  "Static"  and  "Dynamic"  Systems 

Function  space  terminology  provides  a  basis  for  a  precise  definition 
of  the  terms  "static"  and  "dynamic"  as  applied  to  physical  systems.  A  static 
system  is  one  in  which  the  output  is  a  function  of  the  present  input  only. 

Thus,  to  represent  a  general  static  transformation,  it  is  sufficient  to 
expand  the  output  variable  in  terms  of  orthogonal  functions  whose  argument 
is  simply  the  present  value  of  the  input;  it  is  not  necessary  to  utilize  a 
Hilbert  space  description  for  the  input.  A  static  system,  therefore,  is  one 
which  effects  a  transformation  from  a  real  line  to  a  real  line.  In  contrast 
to  static  systems,  in  a  dynamic  system  the  output  depends  not  only  upon 
present  values  of  the  input  but  also  upon  past  values.  Consequently,  a 
dynamic  system  performs  a  transformation  from  a  Hilbert  space  to  a  line. 

Such  systems  are  sometimes  said  to  have  "memory"  or  "energy  storage"  while 
static  systems  are  often  described  as  being  "memoryless".  By  its  nature, 
the  determination  of  the  output  of  a  dynamic  system  requires  a  double 
orthogonal  expansion.  First,  the  past  of  the  input  must  be  expanded  and 

^  Zadeh  also  points  out  that  infinitely  iterated  integrals  can  be  used  with 
appropriate  weighting  functions  to  obtain  an  alternate  representation  for 
nonlinear  operators  which  does  not  involve  expansion  of  the  input. 
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then  the  resulting  coefficients  used  as  independent  variables  for  a  multi¬ 
variable  expansion  representing  the  transformation  to  the  output  variable. 

In  the  remainder  of  this  report,  the  adjectives  "static"  and  "dynamic"  will 
always  be  used  in  accordance  with  the  above  definitions. 

3.3  Experimental  Evaluation  of  the  Characteristics  of  Static  Nonlinear 
Devices 

While  this  report  is  concerned  fundamentally  with  the  characterization 
and  identification  of  nonlinear  dynamic  systems,  some  interesting  simplifica¬ 
tion  of  function  space  methods  occur  when  the  system  under  test  is  either 
static  or  linear.  In  the  case  of  static  nonlinear  systems,  the  output  depends 
only  on  the  present  value  of  the  input  so  an  expansion  of  the  input  in 
orthogonal  functions  is  not  necessary.  The  output  may  be  represented  simply 
as 

co 

y(t)  =  fCx(t)]  *  ^  ai  #i  [x(t)]  (29) 

i=l 

which  is  a  single  rather  than  a  double  expansion.  L.  A.  Zadeh  has  described 
several  suitable  orthonormal  sets  of  functions,  [J2^(x)}  ,  in  a  paper  dealing 
with  static  nonlinearities  which  are  completely  defined  by  their  describing 
functions  (2 )  .  H.  J.  Lory,  et  al.  have  discussed  the  application  of  harmonic 
analysis  utilizing  growing  real  exponentials  to  obtain  the  coefficients  of  a 
Taylor  series  expansion  of  f(x)  (40). 

3.4  Simplification  of  Function  Space  Representations  fox’  Linear  Time 
Invariant  Dynamic  Systems 

When  the  system  under  test  is  linear  and  time  invariant,  only  an  input 
function  space  is  required.  The  transformation  from  the  input  function  space 
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to  the  output  space  need  not  be  expanded  in  orthonormal  functions  since  the 
linearity  of  the  system  demands  that  the  output  be  expressable  as  a  simple 
linear  function  of  the  input  Hilbert  space  coordinates.  For  example, 

Y.  W.  Lee  shows  that 

\ 

GD 

h(t)  -  £  ci  Ai(t)  (30) 

i*l 

where  A^(.t)  are  the  coefficients  of  an  expansion  of  the  part  of  the  system 
input  in  terms  of  orthonormal  Laguerre  functions  and  the  are  coefficients 

characterizing  the  system  (41).  The  orthonormal  Laguerre  functions  used  by 
Lee  are  obtained  from  an  orthogonalization  of  the  family  of  functions 

f  *  (at)n  eat  ,  n  =*  0,1.2,  •••  od  (31) 

n 

Another  example  of  linear  system  representation  utilizing  a  single  expansion 
is  provided  by  T.  P.  Goodman  and  J.  B.  Reswick  (29).  As  previously  mentioned, 
Goodman  and  Reswick  used  a  finite  term  approximation  to  the  expression 

lim  "  N 

y(t)  =  T  -»0  y  T-f  x(t-nT)  h(nT)}  (32) 

NT  =  t  n=0 

In  their  mechanization,  the  values  of  x(t-nT)  are  obtained  from  a  tapped 
delay  line  while  T*h(nT)  represents  potentiometer  settings  weighting  these 
delayed  values.  Other  methods  for  describing  linear  systems  in  terms  of  a 
single  expansion  are  given  in  Lee  (41)  and  T.  Kitamori  (42). 
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3.5  The  Wiener  Theory  of  Nonlinear  Systems 

While  the  idea  of  expanding  operators  over  Hilbert  spaces  is  not  new, 
it  appears  that  only  recently  has  an  explicit  experimental  technique  been 
proposed  for  the  evaluation  of  the  coefficients  in  such  expansions.  A 
specific  method  was  described  by  Norbert  Wiener  in  a  summer  lecture  series 
in  1953-54.  However,  the  first  documentation  of  Wiener's  approach  apparently 
occurred  with  the  publication  of  A.  Bose's  dissertation  in  1956  (3)«  Since 
it  is  felt  that  the  experimental  approach  taken  by  Wiener  represents  the 
only  reasonably  well  explored  alternative  to  the  methods  proposed  in  this 
report,  the  elements  of  the  Wiener  theory  are  summarized  fairly  completely 
in  the  following  paragraphs. 

As  proposed  by  Wiener,  the  input  to  a  system  is  expanded  in  Laguerre 
functions.  These  functions  are  derived  from  the  conventional  Laguerre  poly¬ 
nomials  by  including  the  square  root  of  the  Laguerre  weighting  function, 
e-t,  as  a  multiplying  factor  on  each  polynomial.  Thus,  since  the  nth  Laguerre 
polynomial  is  given  by 


,  ,  ,  (n-1) 

Lnlt)  -TZU’.  e  ftT^nO 


(n-1)  -t"\  tot 
e  j  ,  n=l,2,3 


00 


(33) 


the  n-th  Laguerre  function  may  be  written 


hn(t) 


L  ■  ( t ) 
n 


t  >  0 


t  <  0 


(34) 


The  functions  ^(t)  are  orthonormal  over  [0,ao]with  unit  weighting  function. 


Using  these  functions,  the  past  of  the  input  may  be  expanded  at  any  instant 
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t  >  0 


(35) 


x(-t)  -Yu  h  (t) 
L  n  n 

n-1. 


where 


f*® 

u  -  x(-t)  h  (t) 

n  J  n 


dT 


(36) 


Now  this  equation  has  the  form  of  a  convolution  integral.  Furthermore,  the 
functions,  h^Ct),  are  of  the  same  form  as  the  impulse  response  of  a  linear 
lumped  constant  network  with  a  pole  of  order  n  .  It  is  to  be  expected, 
therefore,  that  the  desired  Laguerre  coefficients  can  be  obtained  continuously 
in  time  by  feeding  x(t)  into  a  bank  of  appropriate  linear  filters.  Wiener 
pointed  out  that  this  is  indeed  the  case.  The  appropriate  filter  transfer 
function  may  be  determined  by  application  of  the  Laplace  transform  to 
equation  33  yielding  (41) 


1  n-1 

s  —  T- 


>.<•>-  (-4) 


(37) 


S  +  2  S  +  2 


This  result  is  illustrated  in  Figure  2.  It  is  interesting  to  note  that  the 

filter  depicted  is  nothing  more  or  less  than  a  low  pass  filter  followed  by  a 

5 

lumped  constant  approximation  to  a  tapped  delay  line  .  This  type  of  filter 
is  easily  constructed  using  standard  analog  computer  elements. 


Since  the  Laguerre  functions  form  a  complete  basis  for  bounded,  continuous, 
square  integrable  functions,  any  nonlinear  system  of  the  class  treated  by 
Wiener  can  be  represented  (for  such  inputs)  by  a  bank  of  linear  filters 
(as  in  Figure  2)  followed  by  a  "zero  memory"  nonlinear  function  generator 
having  the  u^(t)  as  inputs  and  y(t)  as  an  output. 
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In  order  to  fully  exploit  the  possibilities  of  a  Laguerre  function 

expansion  of  the  input  signal,  the  Wiener  theory  requires  that  the  input 

used  to  probe  the  system  response  be  a  broadband  Gaussian  process  such  as 

shot  noise.  With  this  choice,  it  turns  out  that  the  Laguerre  coefficients 

are  themselves  uncorrelated  Gaussian  random  processes  with  equal  variances. 

This  being  the  case,  it  seems  natural  to  expand  the  system  operator  in 

Hermite  functions  since  the  Hermite  polynomials  are  orthonormal  over 

[-CD,  oo]  with  a  Gaussian  weighting  function.  If  T^Cx)  is  the  n-th 
,  _x2 

Hermite  polynomial  (orthonormal  with  e  weighting  function),  then  the 
n-th  Hermite  function  as  defined  by  Wiener  is,  given  by 

2 

x 

<f>n(x)  *  e  2  ^(x)  (38) 


These  functions  are  orthonormal  in  [-oo,  od]  with  unit  weighting  function. 
Weiner  has  shown  that  the  transformation  from  the  Laguerre  coefficient  input 
space  to  the  system  output  can  be  written  as  an  expansion  in  terms  of  Hermite 
functions.  Specifically  (3) 


y(t)  = 


lim 


CD  CD  00 

-CD  l  l'"laij...h  W  +j(U2)  W 

i=l  j=l  h=l 


(39) 


The  coefficients  in  this  expansion,  can  136  determined  by  multiplying 

both  sides  of  the  equation  by  the  appropriate  products  of  (^(x)  functions 
and  averaging  over  [-oo,  od]  .  However,  due  to  Wiener's  judicious  choice  of 


input  signals  and  expansions,  the  necessary  averaging  can  be  accomplished  by 


simply  cross-correlating  the  system  output  with  the  proper  products  of 
Hermite  polynomials.  Thus  (3*38 ) 


s/2  li®  T 

'  <2'>  T  ^oo  h  J.T  »(t)  VV  ty“2>  •••  w  dt  (lM) 

This  can  be  written  more  compactly  by  denoting  the  set  of  subscripts  by 
a  ,  the  product  of  Hermite  polynomials  by  V^Cu),  and  the  time  averaging 
by  the  conventional  bar  symbol  of  statistics.  Thus 

a  =  ( 2n)s/ 2  y(t)  V  (u)  (41) 

a  a 

This  equation  summarizes  the  experimental  part  of  the  Wiener  theory. 

Figure  3  is  a  schematic  representation  of  the  equipment  required  to  carry 
out  the  operations  indicated  (3). 

3.6  Difficulties  Associated  with  Implementation  of  the  Wiener  Theory 

Despite  the  apparent  simplicity  of  equation  4l,  about  all  that  can 
be  said  in  its  favor  is  that  only  a  countable  infinity  of  limiting 
operations  is  involved  in  the  evaluation  of  the  required  coefficients. 

This  is  of  scant  comfort  to  an  investigator  faced  with  the  task  of  actually 
determining  the  characteristics  of  a  real  physical  system.  In  order  to 
actually  make  use  of  the  Wiener  theory,  it  is  necessary  to  truncate  all 
of  these  limiting  operations  both  with  regard  to  measurement  time  and  the 
number  of  terms  taken  in  the  series  for  y(t)  .  So  far  as  is  known  to  the 
writer,  there  has  been  no  analysis  of  the  errors  of  such  truncations. 
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FIG.  3  EXPERIMENTAL  ARRANGEMENT  REQUIRED  FOR  THE  EVALUATION  OF 
WIENER  COEFFICIENTS. 


Such  analysis  would  certainly  be  very  difficult  since  nested  limiting 
operations  of  high  dimensionality  are  involved  in  the  Wiener  theory.  In 
addition  to  this  analytic  difficulty,  there  are  serious  problems  involved 
in  attaining  the  computing  speeds  needed  to  implement  the  experimental 
arrangement  shown  in  Figure  J>.  It  appears  that  a  hybrid  computer  possess¬ 
ing  the  best  features  of  both  analog  and  digital  devices  will  be  required. 
The  extensive  literature  search  undertaken  in  conjunction  with  the  writing 
of  this  report  has  produced  no  reference  relating  an  actual  experiment  of 
the  type  indicated  by  Figure  3. 

The  difficulties  discussed  above  relate  to  dimensionality  and  time 
of  observation.  Even  if  these  problems  should  be  resolved,  there  are  other 
fairly  serious  weaknesses  inherent  in  the  Wiener  theory.  First  of  all, 
there  remains  the  problem  of  obtaining  some  relationship  between  the 
Wiener  coefficients  and  meaningful  system  parameters.  That  is,  from  an 
engineering  point  of  view,  it  would  be  very  desirable  (and  in  many  cir- 
cumstances  essential)  to  invert  the  Wiener  coefficients  to  obtain  the 
parameters  of  the  system  differential  equation.  There  is  no  evidence 
that  this  can  be  done.  Secondly,  the  Wiener  theory  provides  a  basis  only 
for  a  "bench  test"  type  of  experimental  system  identification  somewhat 
analogous  to  frequency  testing  for  linear  systems.  The  experimental 
technique  requires  that  the  system  under  test  be  disconnected  from  its 
normal  operating  signals  and  subjected  to  a  special  test  signal  over  a 
long  period  of  time  and  under  rather  ideal  conditions  of  observation. 

This  is  not  possible  in  many  identification  problems  of  practical  impor¬ 
tance.  Finally,  in  addition  to  all  the  other  stumbling  blocks,  there 
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remains  the  fact  that  the  Wiener  theory  is  not  completely  general.  It 
excludes  the  very  important  classes  of  time  varying  and  unstable  systems 
from  consideration . 

Regardless  of  these  shortcomings,  the  Wiener  theory  appears  to 
stand  as  the  only  general  experimental  technique  for  the  identification 
of  nonlinear  systems  which  has  been  proposed  up  to  the  present  time. 

Whether  or  not  the  theory  is  ever  actually  used  in  its  present  form,  it 
serves  a  very  useful  purpose  in  providing  a  new  conceptual  basis  for  the 
experimental  aspects  of  nonlinear  theory.  The  practical  problems  associated 
with  the  theory  are  certainly  worthy  of  further  attention  . 

4.  PARAMETER  SPACE  METHODS 


4.1  Introduction 

The  Wiener  theory  of  nonlinear  systems  is  applicable  in  situations 
where  a  complete  state  of  ignorance  exists  concerning  the  nature  of  the 
nonlinear  system  under  test.  This  is  rarely  the  case  in  practice.  As  a 
general  rule,  systems  or  devices  subjected  to  experimental  tests  are 
governed  by  physical  principles  which  are  reasonably  well  understood. 

The  uncertainty  in  such  tests  usually  relates  to  the  magnitude  of  various 
effects  rather  than  to  the  basic  mechanisms  operating  to  produce  the 
observed  data.  Even  when  this  is  not  the  case,  the  investigator  is  at 
least  able  to  suggest  several  competing  theories  to  explain  the  observed 
phenomena.  Under  these  circumstances,  it  is  possible  to  construct  a 
finite  parameter  model  for  the  system  under  test  rather  than  an  infinite 
parameter  model  as  proposed  by  Wiener.  This'  approach  will  be  called  a 
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"parameter  space  method"  in  this  report  to  distinguish  it  from  the 
function  space  methods  used  by  Wiener  and  others.  Parameter  space  may  be 
viewed  as  a  generalization  of  the  familiar  "phase"  or  "state"  space  employed 
in  mechanics  for  the  description  of  the  state  of  a  physical  system. 

When  the  object  under  test  is  a  dynamic  system,  a  very  natural 
finite  parameter  characterization  may  be  obtained  by  utilizing  the  system 
differential  equation.  The  previously  mentioned  damped  pendulum  equation 
provides  a  simple  example.  Referring  to  equation  1,  the  pendulum  angle, 

0,  is  governed  by 

10  +  B©  +  Mg i  sin  0*0  (42) 

which  can  be  normalized  to 

c2©  +  c^O  +  sin  0*0  (43) 

The  determination  of  and  c^  along  with  two  initial  conditions 

provides  a  complete  characterization  of  the  system.  This  report  is  bas¬ 
ically  concerned  with  the  formulation  of  an  approach  to  permit  the  inference 
of  parameters  of  this  type  from  unreliable  records  of  the  input  and  output 
of  a  system  under  test. 

As  mentioned  at  the  beginning  of  this  report,  the  idea  of  representing 
a  system  by  a  differential  equation  is  scarcely  a  new  concept.  Indeed  it 
would  seem  to  be  somewhat  paradoxical  to  advocate  a  return  to  differential 
equation  models  in  light  of  the  remarks  made  previously.  The  explanation 
for  this  apparent  regression  lies  in  the  emergence  of  electronic  computers 
as  a  revolutionary  force  in  scientific  and  engineering  methodology. 
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While  it  is  indeed  difficult  to  solve  a  nonlinear  differential  equation 
such  as  equation  43  by  analysis,  it  is  quite  simple  to  obtain  a  solution 
by  electronic  computation.  In  fact,  electronic  analog  computers  are  specif- 
ically  constructed  for  the  purpose  of  solving  high  order  sets  of  nonlinear 
differential  equations  and  are  poorly  suited  to  any  other  task.  Graphical 
and  various  approximate  methods  which  have  been  found  extremely  valuable 
by  human  investigators  are  of  little  or  no  value  to  a  computer.  In 
recognition  of  this  fact,  the  parameter  space  methods  to  be  described  in 
the  remainder  of  this  report  are  entirely  computer  based.  Both  the 
determination  of  the  system  parameters  and  the  evaluation  of  the  resulting 
response  will  be  accomplished  by  completely  automatic  methods. 


4.2  An  Abstract  Comparison  of  Function  Space  and  Parameter  Space 
Characterizations 

In  the  Wiener  theory,  the  input  to  a  system  is  regarded  as  a  trajec¬ 
tory  in  a  Hilbert  space.  Time  appears  parametrically  along  this  curve. 

The  output  at  a  particular  instant  is  obtained  by  the  application  of  a 
transformation  from  the  appropriate  point  in  the  Hilbert  space  to  the 


real  line.  This  transformation  is  completely  described  by  the  infinite 
set  of  Wiener  coefficients,  •  Figure  4  shows  this  situation 

graphically.  In  this  figure,  T^j  stands  for  the  transformation  given 
by  equation  39.  The  distinctive  feature  of  the  expansions  chosen  by  Wiener 
is  that  when  the  trajectory  in  the  input  space  is  produced  by  a  wide-band 


Gaussian  random  process,  the  transformation  from  the  input  space  to  the 
output  space  is  uniquely  determined  and  simply  (conceptually)  computed 
from  simultaneous  observation  of  the  input  and  output  trajectories. 
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In  contrast  to  the  Wiener  theory,  in  the  parameter  space  method 
uncertainties  relating  to  both  the  input  and  the  system  (including  initial 
conditions)  are  described  by  a  finite  number  of  parameters.  Moreover,  due 
to  the  uniqueness  properties  of  solutions  to  specified  differential  equa¬ 
tions,  the  output  variable  is  completely  determined  for  all  time  from  the 
initial  point  in  parameter  space.  Thus,  solution  of  the  differential 
equation  is  equivalent  to  performing  a  transformation,  Tq  ,  from  a  point, 
c  ,  in  a  finite  dimentional  vector  space  to  a  point  in  function  space. 

Figure  5  illustrates  this  relationship  Because  the  transformation, 

Tq(c),  operates  in  a  space  of  finite  dimensionality,  deducing  the  system 
description  is  no  longer  a  problem  in  functional  analysis.  Rather,  in  the 
formulation  to  be  employed  in  this  research  program,  the  system  differential 
equation  is  found  by  utilizing  the  techniques  of  nonlinear  programming. 

This  approach  will  be  explained  in  detail  in  subsequent  reports. 

4.3  Choice  of  a  Metric  for  the  Output  Function  Space 

In  order  to  permit  iterative  techniques  to  be  employed  in  the 
estimation  of  parameter  vectors,  it  is  necessary  to  define  a  distance 
function  or  metric  which  measures  the  distance  between  two  functions. 

For  the  purpose  of  this  investigation,  the  metric  chosen  is  the  conventional 
"Euclidian"  or  L.,  metric 

p2(ylt  y2)  -  f  Cy-^t)  -  y2.(t-)]2  dt  (44) 

J  a 

- —  - -  -  -  -  ■■  —  ■ 

^  It  is  possible  that  the  input  signal,  x(t),  may  be  a  random  process. 

In  that  event,  for  identification  processes  taking  place  in  real  time, 
only  past  values  of  x(t)  will  be  known  so  that  the  output,  y(t),  can  be 
determined  only  up  to  the  present  time.  In  such  circumstances,  the  output 
can  be  thought  of  as  a  trajectory  in  a  space  of  coordinates  associated 
with  a  semi-infinite  "backward  looking"  expansion  rather  than  as  a  single 
point  derived  from  an  expansion  over  [-oo,  +a>]  or  [0,  +a>]  . 
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FIG. 5  PARAMETER  SPACE  DESCRIPTION  OF  A  NONLINEAR 
SYSTEM . 
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With  this  metric,  a  complete  linear  function  space  is  commonly  called  a 
"Hilbert  space".  Note  that  while  p  is  a  functional  in  the  output  space, 
it  is  an  ordinary  function  of  the  parameter  space  coordinates.  This  is  of 
the  utmost  importance  so  far  as  the  computational  aspects  of  system 
identification  are  concerned. 

One  reason  for  choosing  a  Euclidian  metric  is  that  a  considerable 
body  of  knowledge  exists  relating  to  Hilbert  spaces.  However,  an  even 
more  important  reason  from  the  point  of  view  of  this  research  program  is 
that  the  metric  yields  a  set  of  linear  simultaneous  equations  in  the 

process  of  iteratively  deducing  the  system  parameters  from  response  data. 

It  will  be  shown  in  subsequent  reports  that  this  is  the  only  metric  which 
has  this  property.  Since  nonlinear  simultaneous  equations  are  very  difficult 
to  solve  even  by  computer  methods,  this  feature  of  the  L ^  metric  makes  its 
use  almost  mandatory  in  some  of  the  computational  procedures  to  be  described 
in  the  sequel  to  this  report. 


4.4  Experimental  Evaluation  of  Parameter  Space  Coordinates 

The  parameter  space  method  begins  with  the  specification  of  a  finite 
dimensional  parametric  model  for  the  system  under  test.  For  example,  the 
behavior  of  the  system  described  by  equation  43  is  completely  determined 
by  the  specification  of  a  four  dimensional  parameter  vector,  c  : 


(45) 
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In  the  event  that  a  forcing  function  with  certain  unknown  parameters  were 
present,  the  dimensionality  of  this  vector  would  be  increased.  Having 
established  a  model,  experimental  observation  of  the  system  output  may 
begin.  When  a  sufficient  length  of  record  has  been  obtained,  say 
0  <  t  <  T  ,  the  identification  or  inference  process  described  in  the 
following  paragraphs  may  be  initiated. 

For  the  most  general  sort  of  nonlinear  system,  the  parameter 
space  identification  can  be  accomplished  by  an  iteration  procedure 
utilizing  a  metric  in  the  output  function  space.  The  process  begins 
with  an  initial  guess  for  each  of  the  coordinates  of  the  parameter 
vector,  c  .  This  guess  may  be  generated  either  by  a  computer  as  part 
of  the  inference  process  or  it  may  represent  the  best  estimate  of  a 
human  investigator.  Let  this  vector  be  denoted  by  c^.  Furthermore, 
let  the  true  parameter  vector  be  denoted  cq  and  the  observed  system 
response  be  yQ(t)  .  Associated  with  the  vector  c^  will  be  another 
response,  y^(t),  which  can  be  determined  by  a  computer.  When  this  has 
been  accomplished,  the  distance  between  the  two  functions  can  be  computed 
by  evaluating  the  integral: 

p2(y0.  y^  *  J  Cy0(t)  -  y-^t)]2  dt  =  0(^5  yQ)  (46) 

O 

This  notation  emphasizes  the  fact  that  while  p  is  a  functional  in  the 
output  space,  it  can  be  equated  to  a  simple  function,  0  ,  defined  over 
the  parameter„space  for  a  given  set  of  data,  yQ  .  Such  a  function  is 
usually  called  a  "criterion"  or  "objective"  function  in  the  terminology 
of  mathematical  programming  (43,44). 
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Since  0  can  be  evaluated  for  every  c  ,  partial  derivatives  can 
be  determined  in  the  parameter  apace  yielding  a  gradient: 

/  Sc, 


Sc 


M. 


(47) 


M. 

&c 


c  *  c. 


This  gradient  can  serve  as  a  guide  in  choosing  a  new  set  of  parameters 

c  ^  =  C-,  +  Ac  (48) 


such  that 

0(c2  5  yQ)  <  0(cx  ;  yo)  (49) 

Since  for  an  arbitrary  parameter  vector,  c 

0(c  ;  yQ)  =  p2(y,  yQ)  >  0  (50) 

it  follows  that  when  equation  49  is  satisfied  at  every  stage  of  an  iteration, 
the  resulting  sequence,  [0^,  0^,  **♦  0 ^  •••}  must  converge  to  some  limiting 
value,  say  .  Under  ordinary  circumstances,  the  corresponding  sequence 
in  parameter  space  will  also  converge  to  a  limiting  value,  say  c  *  c  • 

When  this  occurs,  the  limit  vector,  cg  ,  provides  an  "estimate"  of  the 
true  parameter  vector,  cq  .  Specific  algorithms  for  obtaining  the  conver¬ 
gence  specified  by  equation  49  are  being  developed  as  part  of  the  continuation 
of  this  research  program. 
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Figure  6  summarizes  the  steps  involved  in  the  experimental  determina¬ 
tion  of  the  parameter  space  coordinates  of  an  unknown  physical  system.  In 
the  actual  implementation  of  such  an  experiment)  all  operations  could  be 
performed  by  a  general  purpose  digital  computer  with  analog  inputs.  However, 
in  many  circumstances  it  would  be  preferable  to  construct  the  block  labeled 
"computer  model"  on  an  analog  computer  under  the  control  of  the  digital 
computer. 

4.5  Local  Minima  and  Non-uniqueness 

The  conceptual  and  experimental  simplicity  of  parameter  space  methods 

is  not  achieved  without  penalty.  Aside  from  the  fact  that  the  method  requires 

considerable  apriori  knowledge  concerning  the  system  under  test,  there  are 

certain  computational  and  mathematical  difficulties  involved.  First  of  all, 

equations  49  and  50  do  not  guarantee  that  ;  yQ)  =  0  .  It  may  turn  out 

that  c  *  c  represents  merely  a  stationary  point  for  the  function  0  in 
e 

the  parameter  space.  That  is,  the  computed  response,  ye(t),  raay  represent 
the  data  very  poorly  even  though  a  small  change  in  any  parameter  makes  the 
fit  even  poorer.  When  this  is  the  case,  correct  identification  of  the  system 
under  test  requires  that  a  more  sophisticated  search  be  carried  out  in  the 
parameter  space  to  locate  the  absolute  minimum  of  0(c ^  ;  yQ)  rather  than 
a  simple  stationary  point.  This  problem  will  also  be  explored  as  part  of 
this  research  program. 

In  addition  to  the  local  minimum  problem,  there  is  a  more  fundamental 
difficulty  associated  with  parameter  space  methods.  Since  It  is  not  assumed 
that  a  special  test  signal  can  be  applied  to  the  system  to  be  identified, 
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FIG. 6  EXPERIMENTAL  DETERMINATION  OF  THE  PARAMETERS  OF  A  DYNAMIC  SYSTEM. 


there  is  no  guarantee  that  the  mapping  shown  in  Figure  5  in  one-to-one* 

It  nay,  in  fact,  be  aany-to-one.  A  single  exaaple  suffices  to  show  this. 
Suppose  that  the  system  under  test  is  an  unforced  linear  time  invariant 
system  with  an  observed  time  response  resulting  entirely  from  unknown 
initial  conditions.  In  such  a  situation,  it  is  possible  to  select  initial 
conditions  which  excite  only  one  normal  mode  of  the  system.  This  being  the 
case,  all  systems  possessing  this  particular  mode  are  indistinguishable 
from  one  another  when  so  excited. 

Abstractly,  the  non-uniqueness  property  of  parameter  space  methods 
may  be  said  to  result  from  the  fact  that  the  metric  used  for  iteration  is 
in  the  wrong  space.  Let  the  ordinary  Euclidian  distance  associated  with  the 
parameter  vector  space  be  denoted  by 

n 

Pp  ("o  ’  1  (coi  '  cki)2  (51) 

i=l 

Then  what  is  really  desired  is  a  computational  algorithm  which  will  guarantee 

pp  (°o  *  °k+l)  <  pp  ('o  •  ck}  (52) 

Unfortunately,  p  (c  ,  c,  )  cannot  be  computed  from  an  observation  of  the 
P  O  K 

response  functions  y  (t)  and  y,  (t)  associated  with  c  and  c,  respec- 

0  K  O  K 

tively.  Instead,  it  is  necessary  to  base  iteration  upon  the  computable 
function  space  distance,  pCy^  i  yQ)»  given  by  equation  46.  Now  in  the 
absence  of  measurement  error,  if  equation  43,  for  example,  is  indeed  an 
exact  description  of  the  system  under  test,  it  does  follow  that  P(V  yk}  ■*’° 
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is  a  necessary  condition  for  ^,(co  ,  c^)  — >  0  ;  i.e.,  for  bo  . 

However,  as  the  example  discussed  previously  shows,  p(yQ  ,  y^)  — *  0  is 
not  sufficient  to  ensure  that  «k  — *  cq  .  In  order  to  make  p(yQ  ,  y^)  -Ha  0 
a  sufficient  condition,  it  is  necessary  to  restrict  the  region  of  parameter 
space  to  be  explored  to  a  subspace  S  in  which  the  mapping  is  one-to-one* 

The  determination  of  suitable  subspaces  by  purely  mathematical  techniques 
is  apt  to  be  extremely  difficult  in  the  majority  of  situations.  It  seems 
more  likely  that  such  regions  would  be  determined  by  a  preliminary  computer 
investigation  of  the  properties  of  the  particular  parametric  model  under 
investigation. 

4.6  Effects  of  Noise,  Measurement  Errors  and  Imprecise  Models 

For  any  one  of  a  number  of  reasons,  the  computer  solution  for  y(t) 
may  fail  to  correspond  exactly  to  the  measured  response  data  even  when  the 
correct  values  of  the  system  parameters  (c  ■  cq)  are  used  in  the  calculation. 
Among  the  major  contributors  to  this  situation  will  be  measurement  errors, 
random  noise  internally  generated  by  the  system  under  test,  and  the  use  of 
a  computer  model  which  ignores  some  of  the  more  subtle  effects  influencing 
the  actual  physical  system.  In  such  circumstances,  the  minimum  value  of 
the  distance  function  0(c  ;  y  )  over  the  permitted  subspace  S  will  be 
greater  than  zero  and  the  parameter  vector,  cg  ,  associated  with  the  minimum 
will  represent  a  "least  squares  estimate"  of  the  true  parameter  value,  cq  . 

The  optimal  quality  of  such  estimates  may  be  shown  in  a  varity  of  circumstances 

(34,45). 
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4.7  An  Alternate  Criterion  function  for  the  Evaluation  of  Paraaeter 
Space  Coordinates 

The  distance  function  0(e  ;  y  )  is  computed  by  an  implicit  procedure 
requiring  solution  of  the  assumed  system  differential  equation.  It  is 
possible  to  use  another  distance  which  can  be  defined  explicitly  in  a 
parameter  space.  Suppose  that  the  experimental  observation  of  a  system 
to  be  identified  is  carried  out  in  such  a  way  that  values  are  obtained  for 
all  derivatives  up  to  the  n^1,  the  order  of  the  systems.  Furthermore, 
suppose  that  the  system  description  can  be  written  in  the  form 

F(a  ,  y)  -  0  (53) 

where  y  is  the  conventional  phase  space  vector,  y  *  (y,  y,  *y,  y^R”^) 

and  a  represents  the  m  remaining  parameter  space  coordinates  (to  be 
determined).  Thus 


where  y  is  an  n  dimensional  vector  and  a  is  a  vector  with  dimensionality 
m  .  If  y  can  indeed  be  measured  without  error,  then  equation  53  must  hold 
for  the  observed  data;  i.e. 

F(ao  ,  yo)  -  0  (55) 

This  equation  must  be  satisfied  for  all  values  of  time  in  (0,  T)  .  Utilizing 
this  fact,  an  objective  or  criterion  function,  <J>(a  i  yQ)  may  he  defined  over 
the  a  space  for  a  given  experimental  record: 
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(56) 


♦(a  5  yQ)  »  J  Y*  (a  ,  yQ(t)^  dt  >  0 

%-  o 

In  light  of  equation  53 t 

4>(a0;  yQ)  -  0  (57) 

This  relationship  may  be  used  to  calculate  aQ  by  iterative  methods  in 
exactly  the  same  way  that  the  function  0(c  ;  y  )  is  used. 

In  an  actual  physical  experiment,  it  will  usually  be  extremely 
difficult  to  directly  measure  all  of  the  system  derivatives.  On  the  other 
„  hand,  attempts  to  differentiate  experimental  data  several  times  are  likely 

to  fail  due  to  the  inevitable  presence  of  small  amounts  of  system  noise 
and  measurement  error.  As  is  well  known,  these  effects  are  accentuated  by 
differentiation  and  may  even  lead  to  derivatives  which  are  unbounded  in  a 
mean  square  sense.  Most  likely,  statistical  estimation  of  the  necessary 
derivatives  will  be  required  to  obtain  values  for  substitution  in  equation 
56.  Whatever  method  is  used  to  obtain  derivatives,  in  the  practical  situa¬ 
tion  the  various  sources  of  error  will  prevent  equation  57  from  holding 
true.  Rather,  as  with  the  function  0  ,  there  will  be  an  estimator,  a@  , 
which  possesses  the  property 

”in  *(3  ,  yQ)  -  ♦(£  ,  yQ)  (58) 

a 

The  vector  is  a  least  squares  estimate  of  a Q  in  a  sense  somewhat 

A  .  r4 

different  from  the  estimate  obtained  using  the  criterion  function  0  .  An 
analog  computer  mechanization  of  equation  58  has  been  described  by  Y.  Kaya 
and  S.  Yamamura  in  a  paper  dealing  with  the  identification  of  linear 
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systems  (46).  Kaya  and  Yaaaoura  used  identical  filtering  on  all  derivatives 
to  obtain  a  relation  derived  from  equation  53  which  would  be  satisfied 
exactly  in  the  absence  of  noise.  The  statistical  optimality  of  their  method 
was  not  discussed. 

the  choice  between  $  or  0  as  a  criterion  function  appears  to 
depend  to  a  large  extent  upon  the  relative  difficulty  of  obtaining  a  suffi¬ 
cient  number  of  derivatives  of  the  system  output  as  compared  to  solving 
the  assumed  system  differential  equation  with  a  computer.  Clearly,  attempts 
to  minimize  4>  by  iterative  nonlinear  programming  methods  will  encounter 
the  same  difficulties  as  occur  in  the  minimization  of  0  .  There  is  one 
special  situation  however,  in  which  the  function  <J>  appears  to  offer 
significant  advantages.  When  the  differential  equation  expressed  by 
equation  53  is  linear  in  the  parameters  a  ,  then  according  to  equation 
56,  <J>  is  a  quadratic  form  in  the  parameter  space.  Moreover,  <|>  is  pos¬ 

itive  definite  and  therefore  has  a  unique  stationary  point  which  is  the 

minimum  existing  at  a  *  a  ,  the  least  squares  estimate  of  a  .  As  a 

e  o 

consequence,  the  non-uniqueness  and  local  minimum  problems  associated  with 
the  general  nonlinear  system  identification  problem  do  not  occur  in  this 
case.  The  nonlinear  system  described  by  equation  43  is  linear  in  the 
parameters  c^  and  c ^  and  so  furnishes  an  instance  in  which  the  above 
statements  apply. 

It  is  worth  noting  at  this  point  that  the  conventional  way  of 
employing  equation  53  when  y  Ct)  is  an  analytic  function  (rather  than  an 
experimental  function)  involves  substituting  yQ  and  its  derivatives  into 
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equation  53  at  m  different  values  of  tine.  This  results  in  m  simultaneous 
equations,  generally  nonlinear,  which  can  be  solved  for  the  m  unknown 
parameters.  While  this  approach  has  also  been  suggested  as  an  experimental 
procedure,  it  does  not  appear  to  be  appropriate  due  to  the  difficulties 
associated  with  differentiation  of  experimental  response  data  (47,48,49). 

4.8  Parameter  Tracking  Servomechanisms 

If  the  minimization  of  a  criterion  function  is  undertaken  in  real 
time  concurrently  with  the  unfolding  of  a  real  physical  process,  then  it  is 
possible  to  devise  computational  procedures  which  permit  parameter  tracking 
for  time  variable  processes  (49,50).  This  approach  permits  a  considerably 
simpler  model  to  be  used  for  the  short  time  description  of  complicated 
systems.  In  most  of  the  work  carried  out  thus  far  in  this  connection, 
analytic  methods,  (i.e.,  analog  computer  methods)  rather  than  algebraic 
methods  of  tracking  have  been  used.  As  a  consequence  of  this  somewhat 
unnatural  restriction,  severe  stability  problems  have  been  encountered  in 
attempting  to  carry  out  actual  experiments.  It  has  been  necessary  to  use 
functions  of  error  rather  than  functionals  in  most  cases  to  achieve  stable 
parameter  determining  loops.  For  example,  M.  Margolis  found  in  the  investiga¬ 
tion  of  a  simple  linear  first  order  system  that  the  parameter  determining 
loop  was  necessarily  unstable  when  an  integral  squared  error  criterion  was 
used  (50).  In  contrast  to  this  situation,  it  is  not  difficult  to  devise 
algebraic  methods  for  parameter  estimation  which  are  always  stable.  The 
development  and  experimental  evaluation  of  such  techniques  is  being  pursued 
in  the  continuation  of  this  research  effort. 
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The  parameter  space  methods  described  in  this  report  are  of  course , 
entirely  applicable  to  time  varying  systems.  The  variation  of  parameters 
may  be  accounted  for  either  explicitly  by  increasing  the  dimensionality  of 
the  parameter  space  or  implicitly  by  evaluating  the  selected  criterion 
function  over  a  sliding  time  interval. 

5.  SUMMARY  AND  CONCLUSIONS 

The  parameter  space  approach  is  believed  to  represent  a  new  point  of 
view  regarding  the  general  problem  of  system  identification.  The  method  is 
restricted  only  by  the  requirement  that  a  parametric  model  of  finite  dimen¬ 
sionality  must  be  provided  for  the  system  under  test.  In  contrast  to  the 
Wiener  approach,  the  parameter  space  method  involves  the  determination  of 
transformations  defined  over  a  finite  dimensional  vector  space  rather  than 
over  a  Hilbert  space.  The  computational  consequences  of  this  difference 
are  of  fundamental  importance;  by  adopting  the  parameter  space  point  of  view, 
a  problem  in  functional  analysis  is  reduced  to  a  nonlinear  programming  problem. 
In  addition  to  the  computational  advantages  gained  by  a  parameter  space 
approach,  there  is  the  further  advantage  that  no  special  test  signals  are 
required  as  in  the  Wiener  theory.  The  model  for  a  nonlinear  system  can  be 
inferred  from  normal  operating  records. 

The  current  phase  of  this  research  is  directed  toward  the  development 
of  explicit  computational  procedures  for  carrying  out  the  proposed  minimiza¬ 
tion  of  error  functionals  by  parameter  space  methods.  This  work  will  provide 
the  basis  for  future  reports. 
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PART  II 


RESEARCH  IN  OPTICAL  COHERENCE 

A  detailed  study  of  certain  aspects  of  the  theory  of  optical 
coherence  has  been  carried  out.  A  major  objective  of  the  study  is  to 
obtain  criteria  for  optimizing  optical  systems  employing  coherent  sources 
and  amplifiers  such  as  lasers.  This  is  the  first  in  a  series  of  studies 
to  investigate  the  effect  of  various  optical  system  devices  such  as 
antennas,  filters,  modulators  and  demodulators,  etc.,  on  the  coherence 
properties  of  optical  signals.  The  present  report  concerns  itself  pri¬ 
marily  with  spatial  and  temporal  filters. 

After  defining  the  weighting  function  and  aperture  field  distribu¬ 
tion  for  a  spatial  frequency  filter,  the  general  transfer  function  for  the 
spatial  filter  is  defined,  with  a  discussion  of  the  approximations  involved. 
The  weighting  function  or  transfer  function  concept  is  then  extended  to 
include  both  spatial  and  temporal  variations.  The  weighting  function  of 
the  system  in  this  sense  is  equal  to  the  response  of  the  filter  to  a 
point  source  in  space  and  a  unit  impulse  in  time.  An  image  equation  is 
then  defined  which  relates  the  image  field  strength  (as  a  function  of 
the  spatial  coordinates  and  time)  to  the  object  field  strength  (as  a 
function  of  the  same  coordinates)  and  the  weighting  function  of  the 
filter,  A  generalized  coherence  function  is  then  defined  as  a  correla¬ 
tion  function  involving  the  spatial  separation  of  two  sources  in  space 
and  the  time  separation  of  two  points  on  a  signal  as  well  as  the  space 
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and  time  variables  themselves.  Ensemble,  time,  and  space  averages  are 
then  defined  and  it  is  shown  that  the  dependence  of  the  coherence  function 
on  the  origin  in  space  and  time  can  be  removed  by  appropriate  time  tyid  apace 
averaging.  This  operation  yields  a  coherence  function  which  is  a  function 
of  the  separation  in  space  and  time  but  not  the  origin  in  these  variables. 

If  the  process  is  stationary  in  these  variables,  no  space  and/or  time 
averaging  is  needed.  A  normalized  correlation  coefficient  or  correlation 
function  can  then  be  defined.  After  these  general  definitions  are  made, 
several  special  cases  of  interest  are  considered.  In  particular,  the 
case  where  the  object  field  strength  is  separable  into  the  product  of  an 
object  field  strength  in  space  and  an  object  field  strength  in  time  is  of 
interest.  The  special  cases  of  combinations  of  perfect  time  coherence  and 
incoherence  and  .perfect  space  coherence  and  incoherence  are  considered  and 
the  coherence  function  evaluated  for  these  cases.  The  image  equation  is 
written  for  each  of  these  cases  and  the  image  function  evaluated.  It  is 
pointed  out  that  the  problem  of  temporal  partial  coherence  or  incoherence 
is  particularly  difficult  to  treat  because  the  spatial  frequency  filter 
transfer  function  is  a  function  of  the  wavelength  and  hence  of  the  tem¬ 
poral  incoherence.  This  means  that  the  transfer  function  of  the  filter 
depends  upon  the  input.  Consequently,  the  general  problem  is  nonlinear 
and  very  difficult  to  treat.  Basically,  what  is  needed  to  handle  this 
problem  are  the  general  solutions  of  Maxwell's  equations,  without  the 
simplifying  assumption  of  a  sinusoidal  time  variation  of  the  driving 
function.  A  first  approximation  to  the  solution  to  this  problem  was 


-57- 


made  by  assuming  that  the  modulation  bandwidth  is  small  compared  to  the 
carrier  frequency,  which  is  almost  always  the  case  in  optics.  This 
assumption  yields  a  quasi-stationary  approach  in  which  the  frequency 
variations  are  sufficiently  slow  so  that  the  response  for  any  given  input 
frequency  will  be  the  steady-state  response  of  the  filtor  ;o  that  frequency. 
The  overall  response  is  found  by  summing  the  responses  to  oach  individual 
frequency. 

This  initial  study  on  the  coherence  properties  of  optical  filters 
will  be  extended  not  only  to  other  devices,  but  will  provide  a  basic 
statistical  description  of  the  processes  to  be  analyzed  in  the  future. 

In  particular,  the  effect  of  coherence  on  information  content  and  rate 
will  be  considered.  Ultimately,  the  optimum  system  will  be  synthesized 
using  the  criterion  of  maximization  of  information  rate  at  the  desired 
output  of  the  system.  The  determination  of  information  rates  will  be 
strongly  affected  by  the  coherence  properties  of  the  signals  involved. 

The  principal  effort  during  this  period  has  concerned  itself  with 
the  effects  of  coherence  on  filtering.  In  particular,  the  effect  of  linear 
spatial  and/or  temporal  frequency  filters  on  optical  waves  has  been  inves¬ 
tigated.  Optical  filters  differ  from  ordinary  low  frequency  filters  in 
that  both  spatial  and  temporal  frequency  effects  must  be  taken  into 
account.  Let  the  weighting  function  G  of  a  spatial  frequency  filter, 
such  as  an  optical  antenna,  be  defined  as  the  response  of  the  filter  to  a 
point  target  in  the  far  field.  It  is  a  function  of  the  image  coordinates 
x'  and  the  y'  ,  the  object  coordinates  x  and  y  ,  and  of  time  t  .  Let 


-58- 


us  temporarily  assume  that  the  time  variation  is  sinusoidal  and  determine 
what  the  spatial  weighting  function  must  be  for  this  case.  The  antenna 
spatial  weighting  function  or  "spread"  function  may  be  related  to  the 

* 

aperture  field  distribution  across  the  spatial  filter.  It  can  be  shown 
that  the  far  field  pattern  arising  from  an  aperture  distribution  F(5,fl) 
with  nearly  uniform  phase  across  the  aperture  is  given  by 

e-jkE 

G(x,y)  =  njk — - —  g(x,y)  (1) 

where 

g(x,y)  *  — i-r  f  F(g,t])  ejk(x5  +  y7])  d|  dT|  (2) 

(2n) 

and  E  *  distance  between  the  aperture  and  a  far-field  point 
k  «  2.n/\ 

The  object  coordinates  x  and  y  are 
x  m  sin  0  cos  0 

y  =  sin  0  sin  0  (2a) 

All  the  object  dimensions  are  given  in  terms  of  spherical  angular 
co  rdinates  whose  origin  is  coincident  with  the  phase  center  of  the 
spatial  frequency  filter.  The  object  field  strength  is  0(x,y,t)  . 

Note  that  the  object  field  strength  in  general  depends  not  only  on  the 
position  of  the  source  in  space  but  also  on  time.  The  time  dependence 

* 

Silver,  S,,  "Microwave  Antenna  Theory  and  Design",  McGraw-Hill,  New 
5'ork,  19^9,  p.  175,  equation  9> 
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nay  rise  from  several  sources.  In  the  first  place,  the  energy  illuminating 
the  object  or  the  radiation  from  a  passive  object  will  in  general  be  time- 
varying  so  that  the  object  field  strength  will  also  be  time-varying.  Or 
if  the  object  is  a  source,  the  source  field  strength  will  be  time-varying. 
In  addition,  the  object  may  be  moving  so  that  it  will  be  located  at  a 
different  position  at  a  later  instant  of  time  yielding  a  time  dependence. 
This  dependence  is  implicit  rather  than  explicit.  The  image  field  strength 
I  depends  upon  the  coordinates  x'  and  y'  along  the  image  surface  as  well 
as  on  time,  where 

x'  =  sin  O'  cos  0' 

y'  =  sin  O'  sin  0'  (2b) 

Note  that  equation  (2)  is  a  Fourier  transform,  that  is,  the  far  field 
pattern  g(x,y)  is  the  Fourier  transform  of  the  aperture  field 
distribution  F(5,Tl)  .  Since  x  and  y  are  spatial  coordinates,  the 
corresponding  transform  variables  k?  and  kTj  have  the  dimension  of  spa¬ 
tial  frequency.  Let  us  therefore  relabel  both  |  and  T]  as 


^  *  0) 

■  2nf 

X 

X 

T]  «  a) 

=  2nf 

1  y 

y 

where  f  and  f  are  spatial  "frequencies"  which  are  the  Fourier 
x  y 

transform  variables  corresponding  to  x  and  y  ,  respectively.  For  the 
special  case  of  nearly  uniform  phase  across  the  aperture,,  the  Fourier 
transform  pair  relating  the  far  field  pattern  G(x,y)  and  the  aperture 
field  distribution  F(?,T|)  is  (except  for  a  constant) 
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g(x,y)  ■  ■  1|-gl  f  F(u  ,  m  )  e^XWx  +  y“y^  du>  dm 

(2k/  ja  x  y  x  y 

JJ  -jk(xu>  +  yu  ) 

F(ux,  a)  )  -  JJ  g(x,y)  e  y  dx  dy 


When  the  simplifying  assumption  of  uniform  phase  cannot  be  made,  the  aperture 

field  distribution  F(u>  ,  qj  )  corresponding  to  the  given  far  field  pattern 

x  y 

must  be  found  by  solving  the  integral  equation 


.  -JKK 

G(x,y)  =  4-  —  F(u  ,  at  )  (cos  0  +  i  *  s) 

H71  K  J  »  x  y  z 


j  k(xu)x  +  yuiy) 


do)  dot 

x  y 


where  i  is  a  unit  vector  normal  to  the  aperture  and  the  s  is  a  unit 
z 

vector  normal  to  the  phase  front.  Since  G(x,y)  is  the  response  of  the 
spatial  filter  to  a  point  source,  the  filter  has  an  impulse  response 
G(x,y)  .  The  transfer  function  of  this  filter  is  then  the  Fourier  trans¬ 
form  of  G(x,y)  or 


T(vV  -  II  $  *7“  [I.  F(vV 


(cos  0  +  I*  •  ^)e 
z 


jk(xux+jcuy) 


i)  dm 

x  yj 


-j(xux+yuy) 


dx  dy 


In  the  special  case  of  nearly  uniform  phase  across  the  aperture,  equation 
(6)  becomes 


Ibid,  p.  173,  equation  8. 


i 


-jwi  +  yuv) 

T(ux,uy)  -  xjk  ?—R—  JJ  g(x,y)  a  y  dx  dy 

-CD 


-jkB 

-  njk  5’(ux,uy)  (7) 


where 


u  *  ko>  (8) 

y  y 

Equation  (7)  shows  that  to  within  a  constant  the  spatial  frequency  transfer 
function  of  the  filter  is  equal  to  the  aperture  field  distribution  over  the 
aperture  of  the  filter,  providing  the  phase  is  reasonably  uniform  over  the 
aperture.  If  the  phase  differs  markedly  from  being  uniform  over  the 
aperture,  then  the  aperture  field  and  transfer  function  are  related  by 
the  more  general  equation  (6). 

In  the  more  general  case  where  the  temporal  as  well  as  the  spatial 
filtering  effect  of  the  filter  is  taken  into  account,  the  overall  weighting 
function  of  the  filter  is  G(x,y,t)  .  The  function  G(x,y,t)  is  equal  to 
the  response  of  the  filter  to  a  point  source  in  space  and  a  unit  impulse 
in  time,  or  in  other  words,  the  response  of  the  filter  is  an  object 
specified  by 


0(x,y,t)  *  6(x,y,t)  =  6(x,y)  6(t) 


(9) 


The  corresponding  spatial  and  temporal  frequency  transfer  function  is 
given  by 
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i 


T(uxtiy») 


(10) 


where 


CD  00 

■  JJ  dx  dy  J  dt  Q(x,y,t) 


J(xu  +  yu  +  «t) 

•  y 


-00  -0D 


u>  ■  2nf  ■  *  kc  (11) 

This  more  general  weighting  function  Q  takes  into  account  both  the  spatial 
and  temporal  filtering  properties  of  the  filter.  In  order  to  find  the 
response  of  the  filter,  or  the  image  field  strength,  from  an  extended 
source  or  object  (designated  by  0(x,y,t)),  the  object  field  strength  0 
is  convolved  with  the  weighting  function,  Q  .  Thus,  when  both  spatial 
and  temporal  frequency  filtering  effects  are  taken  into  account,  the  image 
field  strength  at  the  output  of  a  two  dimensional  linear  additive  spatial 
filter  is  given  by 

1  t 

I(x',y’,t)  =  ff  dx  dy  f  da  0(x,y,a)  G[x'-x,y'-y,t-a]  (12) 

-i 

If  a  3  dimensional  Fourier  transform  is  taken  of  both  sides  of  equation  (12), 
the  transform  relationship  is  given  by 

Jl  [«  ,  U  ,  0)]  m&'iu  ,  U  ,  U>)  T(u  ,  U  ,  Ul)  (13) 

*  J  a  jr  a  j 

where  is  the  Fourier  transform  of  the  image  field  strength  and  is 
the  Fourier  transform  of  the  object  field  strength. 


For  convenience,  let  us  now  restrict  outselves  to  one  spatial 
variable  x(t)  .  In  general,  the  object  field  strength  will  be  random  and 
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will  be  a  function  of  two  independent  variables,  x  and  t  .  The  random 
process  0(x,t)  will  be  completely  defined  by  the  following  2N  dimen¬ 
sional  probability  distribution: 

pCO)  =  p[0(x1,t1),  0(x2,t2),***,0(x1,tN)5  0(x2,t1),0(x2,t2)l»**,0(x2,tN); 
•••••;  OCxj^,^),  0(xN,t2),*-*,  Otx^tjj)]  (14) 

where  0  represents  the  2N  dimensional  row  matrix.  For  convenience 


o<V‘i>  •  °ij 


A  complete  statistical  description  of  the  process  is  usually  obtained  in 
much  fewer  dimensions  than  2N  as  indicated  in  equation  (14),  depending 
on  the  nature  of  the  process.  The  general  "coherence"  function  will  herein 
be  defined  as  the  correlation  function 


R  (Ax,  At,  x,t)  =  0(x,t)  0(x+Ax,  t+At)  =  0(x..,t.)  0(x.,.,t.  ,) 
o  i  j  i+l  j+l 


[r  0(x.  ,t.)  0(x.  .,t.  ..  )  p(0.  .,0.  )  dO.  .  dO.  .  .  - 

JJ  i’  J  l+l’  j+1  r  ij’  l+l, j+1  ij  l+l, j+1 


where 


V* 


xi  *  x 


V  - t  + 


xi+]^  »  x  +  Ax 


The  overscore'  indicates  an  ensemble  average  as  indicated  in  equation  (l6). 
Note  that  the  object  correlation  function  defined  in  equation  (16)  is  a 


I 


function  not  only  of  the  space  and  time  difference  Ax  and  At  respec¬ 
tively,  but  also  of  the  origin  in  space  and  in  time  x  and  t  .  This!  is 
the  most  general  case.  If  Rq  is  explicitly  a  function  of  x  it  is  said 
to  be  space  nonstationary,  if  it  is  explicitly  a  function  of  t  it  is 
said  to  be  time  nonstationary,  if  it  is  explicitly  a  function  of  both  x 
and  t  it  is  said  to  be  both  space  and  time  nonstationary.  If  the  spatial 
bounds  on  the  object  field  strength  0  are  +  X  ,  and  if  the  temporal 
bounds  on  the  object  fields  time  strength  are  +  T  ,  let  us  define  the 
ensemble  average  of  a  function  g  as 

OD 

g[y(x,t),  X,t)  ]  =  J  g[y(x,t),  x,t]  p  [y(x,t) ,  x,t]  dy  (18) 

-CD 

Let  us  define  the  time  average  as 

g[y(x,t),  x,t]  =  ^  ^  I  sb^**^*  x*t^  dt  d^) 

Let  us  define  the  space  average  as 

--------  X 

g[y(x,t) ,  x, t]  =  J  g[y(x,t),  x,t]  dx  (20) 

If  the  process  is  time  stationary,  equations  (l8)  and  (19)  are  equal.  If 
the  process  is  space  stationary  equations  (l8)  and  (20)  are  equal.  If  the 
process  is  both  time  and  space  stationary,  equations  (.18),  (19)  and  (20)  are 
all  equal  to  each  other.  The  space-and-time-averaged  correlation  function, 
or  coherence  function,  is  then 
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H  (Ax, At)  -  R  (Ax,At,x,t)  (21) 

0  0 

If  the  process  Is  both  time-and  s pace -stationary,  then  the  time  and  space 
averaging  indicated  in  equation  (21)  are  superfluous,  since  R  does  not 
depend  explicitly  on  x  and  t  in  this  case.  The  normalized  coherence 
function,  or  correlation  coefficient,  is  defined  as 


Po(Ax,At,x,t) 


Ro<Ax,At,x,t)  -  0(3^., t1:)  0(xi+1,t1+i) 

ao(xiIl*Vl5 


(22) 


where 


0(x,t)  = 

0(x,t)  p  [0(x,t)]  dO 

(25) 

-.op 

,00 

02(x,t)  * 

!  02(x,t)  p  C(x,t)3  do 

(24) 

=  R  (Ax  =  0,  At  =  0,  x,  t) 
o 

oq2  (x,t)  o  0^(x,t)  -  0^(x,t)  (25) 

If  the  process  is  time  stationary  but  not  space  stationary  R^  and  pQ 
are  not  explicitly  dependent  upon  t;  similarly,  they  are  space  stationary 
if  they  are  not  explicitly  a  function  of  x  . 

Let  us  consider  some  special  cases  of  interest.  In  certain  cases, 
the  function  0(x,t)  is  separable  as  follows 

0(x,t)  »  0(x)  0(t) 
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(26) 


The  aost  common  examples  of  situations  satisfying  equation  (26)  are  when 
the  object  radiation  or  reradiation  is  monochromatic  or  the  source  is  a 
point  source,  or  both.  Thus,  for  a  monochromatic  source,  equation  (26) 
becomes 

0(x,t)  «  0(x)  cos  a)Qt  (27) 

where  o)Q  is  the  frequency  of  the  monochromatic  radiation.  In  this  case 
equation  (16)  becomes 


K 


i 


R  (Ax,At,x,t)  =  0(x. )  0(x.  .. )  cos  u>  t.  cos  u>  t..T  (28) 

o  l  i+i  o  j  o  J+l 

Note  that  the  cosine  function  is  not  affected  by  the  ensemble  average  since 
it  is  nonrandom.  Furthermore,  0(x,t)  is  not  stationary  in  time,  although 
it  has  been  assumed  that  the  object  is  stationary  in  the  space  variable. 
Therefore,  a  further  average  must  take  place  in  time  yielding 

R  (Ax, At)  *  Q(x.)  0(x.  ,)  cos  u>  t .  cos  w  t .  .. 
o  l  i+i  °  3  o  j  +x 

=  0(x. )  0(x.  , )  cos  at  At  (29) 

1 _ 1+1  o 

2 

Note  that  the  coherence  function  in  equation  (29)  is  periodic  in  x  . 
Coherence  or  correlation  functions  which  are  periodic  in  t  are  said  to 
be  "perfectly  time  coherent"  since  a  time  function  is  perfectly  predictable 
from  a  knowledge  of  its  value  at  any  given  time.  This  is  true  of  periodic 
functions  in  general,  since  they  may  be  represented  as  a  Fourier  series, 
therefore  their  correlation  functions  are  a  sum  of  cosine  terms  of  the  form 
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occurring  in  equation  (29).  Any  periodic  or  predictable  wave fora  has 
zero  bandwidth  occupancy,  in  the  sense  that  although  there  are  many 
frequencies  present,  the  set  of  spectral  frequencies  is  a  countable  set 
of  aeasure  zero.  In  other  words,  the  spectrua  consists  of  delta  func¬ 
tions  at  the  various  discrete  frequencies  involved  with  a  net  zero  aeasure 
for  the  bandwidth  occupancy.  From  henceforth,  then,  a  time  function  is 
said  to  be  perfectly  time  coherent  if  its  coherence  function  is  a  periodic 
function  in  r  . 

Now  let  us  consider  the  case  where  the  object  is  a  point  target  in 
space.  Equation  (26)  then  becomes 

0(x,t)  «  A6(x-xq)  0(t)  (30) 

where  A  is  a  constant.  For  convenience  assume  the  process  to  be  time 
stationary;  then  the  coherence  function  becomes 


R0(Ax,At,x,t)  «  A26(xj-xq)  6(xi+1-xQ)  0(tj)  0(tj+1>  (31) 

The  portion  of  equation  (31)  involving  space  variables  is  not  random  and 
therefore  a  space  average  is  needed  to  remove  the  space  dependence.  Thus 


_ _  X 

R  (Ax, At)  «  A2  0(t,)0(t,  .)  f  6 ( x-x  ) 6 ( x-x  +Ax )  dx 

O  j  j+1  Xr»CQ  2X  J O  O 


where 


A2  0(t.)0(t,  )  J1® 

j  j+1  X-*  DO 

a  =  KA“  otV°<W 

for  Ax  *  0 

-  0 

for  Ax  /  0 

l'im 

-  K  for  Ax  -  0 

X-GD 

2X 

■  0  for  Ax  t  0 

(32) 


(32a) 


! 


:**  * 


i 


The  llffliting  operation;  in  equation  (32a)  is  indeterminate  in  that  both 
the  numerator  and  the  denominator  approach  infinity.  Assuming  that  in 
the  limit  this  ratio  approaches  a  constant  K  for  Ax  ■  0  ,  the  space 
averaged  coherence  function  is  a  constant  for  Ax  *  0  and  is  zero  for 
Ax  /  0  .  This  is  the  limiting  case  of  a  periodic  function  of  a  space 
variable  where  the  period  approaches  infinity  or  the  repetition  frequency 
approaches  0  .  The  function  is  then  said  to  be  "perfectly  space  coherent". 
In  general,  if  the  coherence  function  is  a  periodic  function  of  the  space 
separation  Ax  ,  it  is  said  to  be  perfectly  coherent  in  that  variable. 

An  example  of  a  case  where  the  process  is  both  space  and  time  coherent  is 

G(x,t)  »  Mix):  cos  u>o t  (33) 

In  this  case,  the  coherence  function  is: 

KA2 

R  (Ax  *  0,  At)  ■  cos  (,)  At  (34) 

o  do 

Note  that  equation  (34)  is  periodic  in  both  the  space  and  the  time  variable 
and  therefore  is  both  space-and  time-coherent. 

We  have  now  considered  a  case  where  the  variables  may  be  perfectly 
space  coherent,  perfectly  time  coherent  or  both.  Let  us  now  consider  the 
opposite  extreme  where  the  functions  may  be  perfectly  space  incoherent, 
perfectly  time  incoherent  or  both.  Let  us  again  assume  that  the  object 
is  separable  in  the  sense  of  equation  (26).  If  the  object  is  perfectly 
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space  incoherent  ,  the  correlation  bat ween  x^  and  x^  is  zero  except 
when  x^  equals  x^  A  similar  argument  holds  for  perfect  time  incoherence 
Thus,  the  coherence  function  fora  space-incoherent  object  is 

R  ’(Ax,At)  *  X^Ax)  OU^jo(t^>  -  (35) 

In  equation  (35),  is  proportional  to  the  spectral  density  of  0(x)  . 

Similarly,  for  a  temporally  incoherent  source,  the  coherence  function 
becomes 

R  (Ax, At)  -  0(x.  )0(x.)  X_6(At)  (36) 

O  1  c.  d. 

where  is  the  temporal  spectral  density.  If  the  source,  or  object,  is 

both  spatially  and  temporally  incoherent,  the  coherence  function  becomes 

R  (Ax, At)  =  X  K_  6(Ax)6(At)  (37) 

O  X  c. 

Thus,  if  the  function  is  space  incoherent  it  will  involve  a  delta  function 
with  argument  Ax  if  it  is  temporally  incoherent,  it  will  involve  a  delta 
function  with  argument  At  and  if  it  is  both  spatially  and  temporally 
incoherent  it  will  involve  the  product  of  the  two  delta  functions.  In  the 
most  general  case  defined  by  equation  (1'6),  the  coherence  function  will 
be  neither  periodic  nor  involve  a  delta  function  in  either  of  the  variables 
under  these  circumstances  the  object  is  said  to  be  partially  coherent. 

It  is  important  to  consider  the  coherence  properties  of  the  image 
of  a  spatial  and/or  temporal  frequency  filter  under  these.-' circumstances . 
First  let  us  consider  the  case  whera  the  object  is  perfectly  space  coherent 
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and  time  coherent.  Then,  substituting  equation  (33)  in  equation  (12) 
and  limiting  the  equation  to  one  space  variable  only, 


I* 

I(x* ,t)  m  [  dx  [  da  0(x,a)  Q(x,-x,t-a) 
J-1  “o 


f  dx  f  da6(x)  cos  u>  a  G(x'-x,t-a) 
J-1  Jo  0 


«  J  G(x',t-a)  cos  u^a  da  (38) 


In  a  device  which  is  intended  to  be  a  spatial  frequency  filter  only,  such 
as  an  antenna,  it  is  common  to  try  to  make  the  device  as  independent  of 
temporal  variations  as  possible..  If  this  is  the  case,  the  weighting 
function  of  G(x,t)  is  approximately  independent  of  t  ;  and  consequently 
equation  (38)  becomes 


I(x',t) 


dx  6(x)  G(x'-x)  cos 


U)  t 
o 


dx  =  G(x' )  cos 


b)  t 
o 


(39) 


The  image  coherence  function  for  equation  (39)  is 


Rj(Ax' ,At,x' ,t)  «  I(x\t)  Kx'+Ax'  ,t+At) 

*  G(x' )  G(x4Ax')cos  m  t  cos  w  (t+At)  (40) 

o  o 

Averaging  with  reBpect  to  both  space  and  time, 


R-(Ax'  ,At)  »  — x'~  % izLttell  cos  o>  At  (41) 

£  do 

Since  the  time-averaged  coherence  function  in  equation  (4l)  is  periodic 
in  At  ;  the  image  is  perfectly  time  coherent.  The 


1 


r'.. 


coherence  function  will ,  in  general ,  be  space  coherent  in  spite  of  the 
bandwidth  limiting  effect  of  the  spatial-frequency  filter,  since  the 
spatial  signal  response  is  known  for  all  x'  and  is  not  random.  In  the 
case  of  the  image  represented  by  equation  (38) ,  the  temporal  response 
will  be  perfectly  time  coherent  but  its  coherence  function  will  not  be 
periodic  with  finite  period  until  the  steady  state  is  reached.  When 
the  transient  terms  have  died  out,  then  the  steady-state  response  of  a 
linear  temporal  filter  to  a  cosine  wave  is  another  cosine  wave  changed  in 
amplitude  and  shifted  in  phase.  The  steady-state  is,  of  course,  temporally 
coherent.  Spatial  coherence  should  not  be  expected  in  the  sense  that  the 
spatial  coherence  function  will  be  periodic  with  finite  period  or  constant 
because  there  is,  in  effect,  no  equivalent  steady-state  phenomenon  in  the 
spatial  domain.  The  reason  for  this  is  that  the  bounds  of  integration  on 
the  spatial  domain  are  from  -1  to  +1  and  the  concept  of  steady-state 
loses  its  significance.  Let  us  assume  the  special  case  where  the  weighting 
function  is  separable  into  the  product  of  a  space  weighting  function  and  a 
time-weighting  function,  or 


G(x,t)  -  G  (x)  G  (t)  (42) 

X  t 

For  a  spatial  filter,  every  effort  is  made  to  make  G^tt)  equal  to  one  and 
likewise  in  a  temporal  frequency  filter  every  effort  is  .made  to  make  G^Cx) 
one.  The  latter  problem  seldom  comes  up  because  temporal  filters  are 


usually  used  on  voltages  which  are  pure  time  functions.  Now  let  us  consider 
the  image  function  for  the  case  of  pure  spatial  coherence  only.  Then, 


The  coherence  function  for  this  case  is 


RI(Ax‘,At,x',t)  -  Qx(x')  ax(x*+Ax') 


t  t  _ _ 

J  1  d6  J  2  dp  CXt^-a)  0(t2-p)  Qt(a)  Gt(p) 


o  o 


.t,  .t 


Gx(x-)  Gx(x'+Ax);  J  1  do  J  2  dp  RQ(At4a-p,t)  Gt(a)  GtU)  (44) 


0  0 


In  general,  the  image  in  this  case  will  be  coherent  in  neither  space  nor 
time.  Now  let  us  consider  the  coherence  properties  of  the  image  when  the 
object  is  coherent  in  time  but  not  in  space. 

!  t 

l(x',t)  =  j  0(x)  G^Cx'-x)  dx  J  cos  uio(t-a)  G^Ca)  da  (43 

-1  o 

The  image  coherence  function  in  this  case  becomes 

1 

Rj(Ax'  ,At,x'  *t)  •  |J  0(x)0(y)  G^tx'-x)  G^x'+Ax'-y)  dx  dy  *■ 


.t,  „t~ 


f  1  da  f  2  dP  cos  u)  (t^-a)  cos  <d  (t_-P)  G  (a)  G  (p) 
,)  J  0  1  o  d  t  t 

O  0 

1 

■  JT  Ro(x-y)  G^Ix'+Ax-y)  G^Cx'-x)  dx  dy  • 


„t,  „t 


|  J  1  da  2  dP[cos  u>o(t1+t2-a-p)  +  cos  wo(At+a-P)]  Gt(a)  Gt(p)  (46) 


o  o 


Thus,  even  though  the  object  is  perfectly  correlated  in  time,  the  image 
coherence  function  will  not  be  periodic  with  finite  period  except  in 
the  steady  state.  Notice  that  in  all  of  these  integrations 


) 


the  bounds  on  ths  space  dimensions  of  the  object  ere  -1  to  1.  The 
corresponding  analogous  bounds  on  t  should  be  zero  to  infinity.  It 
is  interesting  that  this  definition  of  the  possibility  of  a  transient 
coherence  function,  that  is,  one  in  which  coherence  is  measured  after  the 
signal  has  been  applied  for  only  a  finite  period  of  time.  It  is  clear 
that  if  the  object  covers  a  finite  solid  angle  instead  of  the  4x  radians 
required  to  take  into  account  all  of  the  volume  of  space,  the  bounds  on 
the  space  dimension  would  be  less  than  from  -1  to  1. 

Now  let  us  consider  the  opposite  extreme,  namely  pure  space  incoherence, 
pure-time  incoherence  and  both  time  and  space  incoherence  jointly.  Assume 
that  the  object  is  separable  in  the  sense  of  equation  (26)  and  that  the 
weighting  function  is  separable  in  the  sense  Of  equation  (42).  The  image 
field  can  then  be  written  as 

f  1  Pt 

I(x' ,t)  *  dx  0(x)  G  (x'-x)  da  0(a)  G  (t-a)  (47) 

J-1  x  -'0 

If  the  object  field  strength  is  purely  incoherent  in  both  space  and  time, 
then 


0(xx)  0(x2)  -  Kx  6(x2-x1)  =  Kx  6(dx) 

(48) 

0(tx)  0(t2)  «=  K2  6(t2-tx)  -  K2  6 (At)' 

(49) 

The  corresponding  coherence  function  for  I  is  then 
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the  bounds  on  the  space  dimensions  of  the  object  are  -1  to  1.  The 
corresponding  analogous  bounds  on  t  should  be  zero  to  infinity.  It 
is  interesting  that  this  definition  of  the  possibility  of  a  transient 
coherence  function,  that  is,  one  in  which  coherence  is  measured  after  the 
signal  has  been  applied  for  only  a  finite  period  of  time.  It  is  clear 
that  if  the  object  covers  a  finite  solid  angle  instead  of  the  4n  radians 
required  to  take  into  account  all  of  the  volume  of  space,  the  bounds  on 
the  space  dimension  would  be  less  than  from  -1  to  1. 

Now  let  us  consider  the  opposite  extreme,  namely  pure  space  incoherence, 
pure-time  incoherence  and  both  time  and  space  incoherence  jointly.  Assume 
that  the  object  is  separable  in  the  sense  of  equation  (26)  and  that  the 
weighting  function  is  separable  in  the  sense  of  equation  (42).  The  image 
field  can  then  be  written  as 

i  t 

i(x' ,t)  *  I  dx  0(x)  G  (x'-x)  da  0(a)  G  (t-a)  (47) 

J-1  x  Jo  1 

If  the  object  field  strength  is  purely  incoherent  in  both  space  and  time, 
then 

0(X  )  0(x2)  -  Kx  6(x2-Xl)  =  Kj  6 (dx)  (48) 

0(tx)  0(t2)  =  K2  6(t2-t1)  -  K 2  6(4t)  (49) 

The  corresponding  coherence  function  for  I  is  then 
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1 

RI(Ax,.At,x',t)  -  JJ  dx  dy  0(x)0(y)  G^x'-x)  Q^x'-y)  • 

-1 

J  1  da  J  ?  dp  O(t-a)  O(t-p)  Qt(a)  0^(0) 

*  0  *  o 


dx  Kx 


da  K2  Qt2(a) 


(50) 


It  is  interesting  to  consider  the-  case  of  a  spatially- incoherent  Or  diffuse 
source  which  emits  monochromatic  radiation.  In  this  case  the  image  field 
strength  becomes 


r  l  Pt 

I(x',t)  *  dx  0(x)  G  (x'-x)  da  G  (t-a)  cos  u  a  (51) 

x  Jo  0 

If  the  temporal  characteristic  of  the  spatial  frequency  filter  is  a  very 
broad  compared  to  the  spatial  attenuation  characteristic, 


G  (t)'  -  6(t) 


The  image  field  strength  in  this  case  is 


I(x' ,t)  =  cos  uot  j  dx  0(x)  G^x'-x) 


(52) 


and  the  coherence  function  is 


R-^Ax'  tAt,t,x' )  ■cos  u)Qt  cos  u»o(t+At)  JJ  dx  dy  • 

-1 

0(x)0(y )  Gk(x'-x)  Gx(x'+Ax'-y)  - 

r  1 

K,  cos  a)  t  cos  to  (t+At)  G  (x'-x)  G  (x'+Ax'-x) 
1  o  o  J  x  x 


dx 


(53) 
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The  corresponding  mean  square  value  for  the  Image,  sometimes  called  the 
image  intensity,  is 

I^(x' ,  t)  *  K  coa2  «  t  f  G^x'-x)  dx  (54) 

O  X 

Similarly,  when  the  object  field  strength  is  spatially  coherent  and 
temporally  incoherent,  for  example,  when  the  object  is  a  point  source 
emitting  white  noise,  the  image  field  strength,  coherence  function  and 
intensity  or  mean,  square  value  are,  respectively, 


I(x' ,t)  ■  G  Cx')  0(a)  G.(t-a)  da 

X  J  X 


(55) 


RjCAx'.dt.x'.t)  •  Gx(*'>  G  (x'+dx')  I*  1  da  f2  dp 

A  J0  J0 


O^-a)  0(t2-p)  Gt(a)  Gt(p)  = 

Gx(x')  Gx(x'+Ax')  k2  |  Gt(a)  Gt(dt+a)  da 

I2  «  G2  (x')  K-  r  G2  (a)  da 
x  2  J  t 


(56) 


(57) 


Equations  (54)  and  (57)  illustrate  the  often-mentioned  fact  that  image 
intensity  is  proportional  to  the  integral  of  the  square  of  the  weighting 
function  of  the  spatial  filter  for  incoherent  sources. 


Before  the  advent  of  the  laser,  most  optical  problems  made  use  of 
intensity  relationships  such  as  those  in  Equations  ($4.)  and  (57) ;  With  - 
the  possibility  of  nearly  coherent  or  at  least  partially  coherent  light, 
however,  coherent  optics  must  be  considered,  and  their  coherence  properties 
can  be  measured  and  evaluated  as  has  just  been  illustrated. 

It  will  be  observed  that  this  weighting  function  G(x,t)  which  is 
the  transform  of  the  spatial  frequency  filter  transfer  function  in  equation 
(7)  is  a  function  of  the  wavelength  (k  *  2rt/X)  .  Thus,  the  spatial 
weighting  function  is  dependent  on  the  temporal  frequency  of  the  input. 

This  means  that  G  in  equation  (1)  is  a  function  of  X  and,  consequently, 
so  is  the  image  field  strength  I  .  This  frequency  dependence  must  be 
taken  to  account  when  the  input  or  object  field  strength  is  temporally 
partially  coherent,  that  is,  has  a  finite  temporal  bandwidth.  In  such  a 
case,  the  time-varying  portion  of  the  object  field  strength  can  iso  be 
written  as  a  function  of  wavelength  or  frequency  by  means  of  the  Fourier 
integral  theorem.  Thus, 


°t<t) 


i  r  £*B 

2n  -2nB 

i.  r2*B 

2n  -2tiB 


0^(,u>)  e^^  dm 

doi  J  ”  O^Cz)  Z^dz 


where 


u>  • 


2tic 

X 


(58) 


=  kc 


(59) 
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0^(t)  is  then  a  function  of  the  spectral  bandwidth  B  ,  and  if  the 
spectral  bandwidth  consists  only  of  discrete  lines ,  and  the  integral  in 
equation  (58)  reduces  to  a  Fourier  series  which  involves  the  various 
frequency  components  in  the  series.  Strictly  speaking,  the  basic  Fourier 
transform  pair,  equations  (7)  and  (8)  were  derived  from  Maxwell's  equations 
on  the  assumption  that  the  time  variation  of  the  driving  function  in  the 
wave  equation  was  sinusoidal  and  is  only  a  steady  state  solution.  To  be 
strictly  rigorous,  these  equations  should  be  solved  for  the  case  when  the 
input  is  not  sinusoidal,  but  rather  as  represented  in  equation  (58),  that  is, 
when  the  driving  function  for  Maxwell's  equations  is  not  sinusoidally 
time  dependent,  but  rather  a  spectrum  of  frequencies.  The  image  time  func¬ 
tion  in  the  steady  state  can  be  found  as  the  sum  of  the  solutions  to  the 
various  frequency  components  emitted  by  the  object  by  employing  equations 
(7)  and  (8)  for  every  frequency  component  and  then  summing  the  outputs, 
providing  tlfe  filter  is  linear.  This  procedure  does  not  in  general  lead 
to  a  closed  form  solution  and  furthermore  does  not  yield  the  transient 
response.  If  the  modulation  bandwidth  is  very  small  compared  to  the  carrier 
frequency,  which  is  almost  always  the  case  in  optics,  the  weighting  function 
of  a  spatial  frequency  filter  can  be  thought  of  as  being  the  response  to  an 
input  sinusoid  which  is  frequency  variant.  This  is  analogous  to  the  so- 
called  quasi-stationary  approach  used  in  f-m.  It  is  assumed  that  the 
frequency  variations  are  sufficiently  slow  so  that  the  response  for  any 
given  input  frequency  will  be  the  steady  state  response  of  the  filter  to 
that  frequency  from  which  the  over-all  response  is  found  by  summing  the 
responses  to  each  individual  frequency.  This  sum  becomes  an  integral  in 
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the  limit  of  &  continuous  frequency  variation.  With  these  approximations 
in  mind,  employing  equations  (7)  and  (8),  the  response  of  the  filter  to  a 
point  source  (for  example,  at  the  origin  in  space)  and  a  sine  wave  of 
frequency  u>  in  time  is 

H(x,y,t,a>)  ■  H(x,y,u>) 

j  |  (x?+11y) 

■  J2fer  v- 1.  *  ■>«  dT> 

A 

,  .  ,  0)R 

jut  ~i  T 

n.lme  e _ _  g(x,y,u>)  (60) 

“  c  R 

In  order  to  apply  equation  (12),  however,  the  function  G(x,y,t) 
which  is  the  response  to  a  unit  impulse  in  both  time  and  space  must  be 
found.  Since  the  temporal  weighting  function  is  simply  the  integral  of 
the  steady  state  response  to  a  sinusoid  in  time  overfall  frequencies,  the 
function  G  can  be  found  from  H(x,y,t,w)  in  equation  (60)  as  follows: 

G(x,y,t)  ■  T  H(x,y,u)  e^wt  du>  (6l) 

-  CD 

Equation  (12)  then  becomes 
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where 

r  jf  (x?+yff) 

g(x,y,«)  -  F(5,T|)  e  d§  d7] 

4*Z  JA 


(62) 


(63) 


If  the  modulation  bandwidth  is  extremely  small  compared  to  the  carrier 
frequency,  to  a  first  approximation  the  weighting  function  in  equation 
(60)  becomes  independent  of  frequency  u  in  the  sense  that  u>  i6  equal 
approximately  to  a  constant.  In  this  case,  the  spatial  weighting  function 
is  independent  of  time,  that  is,  there  is  no  temporal  filtering  by  such 
a  filter.  The  evaluation  of  the  integral  (62)  is  clearly  difficult. 

Even  equation  (62)  does  not  apply  if  the  modulation  bandwidth  of  the 
object  spectrum  is  an  appreciable  fraction  of  the  carrier  frequency. 

When  these  expressions  are  invalid,  the  general  expression  for  the  image 
is  even  more  complicated.  Because  of  the  difficulties  associated  with 
time  incoherent  signals,  the  principal  effort  for  this  report  has  been 
with  temporally  coherent  sources  but  with  any  degree  of  spatial  coherence. 


Work  is  currently  under  progress  to  evaluate  information  rates  under  these 
conditions.  Once  the  information  rate  has  been  evaluated,  it  can  be 
maximized  with  respect  to  the  variable  parameters  of  the  filter,  thereby 
optimizing  the  filter. 

Work  for  the  forthcoming  period  will  include  an  investigation  of 
the  difficulties  associated  with  temporal  incoherence  and  an  attempt  to 
evaluate  information  rates  for  signals  which  are  partially  time  coherent 
as  well  as  partially  space  coherent.  Investigation  will  continue  into 
the  behavior  of  devices  other  than  filters  in  an  optical  system,  and 
particularly  the  effect  of  these  devices  on  coherence  and  information 
rate.  Included  in  this  study  will  be  the  effect  of  receiver  noise  added 
after  spatial  frequency  filtering,  that  is,  in  addition  to  the  background 
noise  present  in  the  image,  The  effect  of  nonlinear  operations  will  also 
be  considered,  especially  simple  nonlinear  algebraic  operations  such  as 
encountered  in  modulators  and  demodulators. 
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